(1) This submission provides a pair of ASCII files to be "sourced" into S. The resulting S functions, named filter and forecast, implement Kalman filter recursions using the EM algorithm. Reference: "An approach to time series smoothing and forecasting using the EM algorithm." R. H. Shumway and D. S. Stoffer J. Time Series Anal. Vol. 3, No. 4 253-264 (2) The information I would like to share: filter <- function(numobs,y,A,Sigma0,mu0,Phi0,Q0,R0,tol=0.01,max.iter=20,verbose=T){ # y[,t] rx1 # A[[t]] rxp # Sigma0 pxp # mu0 px1 # Phi0 pxp # Q0 pxp # R0 rxr # tol, log likelihood convergence tolerance # EM estimation of Kalman filter # Reference: "An approach to time series smoothing # and forecasting using the EM algorithm." # R. H. Shumway and D. S. Stoffer # J. Time Series Anal. Vol. 3, No. 4 253-264 # y(t) = A(t)%*%x(t) + v(t), for t=1,...,T # x(t) = Phi%*%x(t-1) + w(t), for t=1,...,T # y rx1 observation vector # A rxp design matrix # E[v(t)]=0, R = E[v(t)%*%t(v(t))] rxr observation noise covariance # x px1 state vector # Phi pxp transition matrix # E[w(t)]=0, Q = E[w(t)%*%t(w(t))] pxp state noise covariance # initial vector x0, E[x0] = mu, Sigma = E[(x0-mu)%*%t(x0-mu)] # Notation: # sxt = E[x(t)|y(1),...,y(s)] # sPtu = E[(x(t)-sxt)%*%t(x(u)-sxu)|y(1),...,y(s)] N <- numobs mu <- mu0 Phi <- Phi0 Q <- Q0 R <- R0 iter <- 0 cvg <- 1 while((iter <= max.iter) && (cvg > tol)){ # gain K pxr and covariances sPtt, tPtt sPtt <- vector("list",N) K <- vector("list",N) tPtt <- vector("list",N+1) tPtt[[1]] <- Sigma0 sxt <- vector("list",N) txt <- vector("list",N+1) txt[[1]] <- mu eps <- vector("list",N) Sigma <- vector("list",N) qf <- seq(N)*0 det <- seq(N)*0 # forward recursions for(n in seq(N)){ sPtt[[n]] <- Phi%*%tPtt[[n]]%*%t(Phi) + Q K[[n]] <- sPtt[[n]]%*%t(A[[n]])%*%solve(A[[n]]%*%sPtt[[n]]%*%t(A[[n]])+R) tPtt[[n+1]] <- sPtt[[n]] - K[[n]]%*%A[[n]]%*%sPtt[[n]] # filter estimators sxt, txt sxt[[n]] <- Phi%*%txt[[n]] txt[[n+1]] <- sxt[[n]] + (K[[n]]%*%(y[,n]-(A[[n]]%*%sxt[[n]]))) # innovations eps eps[[n]] <- y[,n] - A[[n]]%*%sxt[[n]] Sgm <- A[[n]]%*%sPtt[[n]]%*%t(A[[n]]) + R Sigma[[n]] <- (t(Sgm)+Sgm)/2 # To ensure its symmetric! qf[n] <- t(eps[[n]])%*%solve(Sigma[[n]])%*%eps[[n]] det[n] <- prod(eigen( Sigma[[n]] )$values) # Det for real-valued matrices! # det[n] <- prod(diag(chol( Sigma[[n]] )))^2 } # backward recursions J <- vector("list",N) Txs <- vector("list",N+1) Txs[[N+1]] <- txt[[N+1]] TPss <- vector("list",N+1) TPss[[N+1]] <- tPtt[[N+1]] st0 <- vector("list",N) ss0 <- vector("list",N) St0 <- matrix(0,nrow=dim(TPss[[N+1]])[1],ncol=dim(TPss[[N+1]])[2]) Ss0 <- matrix(0,nrow=dim(TPss[[N+1]])[1],ncol=dim(TPss[[N+1]])[2]) # smoother estimators Txs for(n in seq(from=N,to=1)){ J[[n]] <- tPtt[[n]]%*%t(Phi)%*%solve(sPtt[[n]]) Txs[[n]] <- txt[[n]] + J[[n]]%*%(Txs[[n+1]]-sxt[[n]]) # covariance TPss TPss[[n]] <- tPtt[[n]] + J[[n]]%*%(TPss[[n+1]]-sPtt[[n]])%*%t(J[[n]]) # Q(,) calculations st0[[n]] <- TPss[[n+1]] + Txs[[n+1]]%*%t(Txs[[n+1]]) ss0[[n]] <- TPss[[n]] + Txs[[n]]%*%t(Txs[[n]]) St0 <- St0 + st0[[n]] Ss0 <- Ss0 + ss0[[n]] } loglik <- -.5*sum(log(det)) - .5*sum(qf) if(-2*loglik < 0){ out <- list(mu, Phi, Q, R) names(out) <- c("mu", "Phi", "Q", "R") warning("-2*loglik < 0") return(out) } if(verbose==T) print(-2*loglik) # EM algorithm # covariance TPsr TPsr <- vector("list",N) p <- dim(Phi0)[1] TPsr[[N]] <- (diag(p) - K[[N]]%*%A[[N]])%*%Phi%*%tPtt[[N]] for(n in seq(from=N,to=2)){ TPsr[[n-1]] <- tPtt[[n]]%*%t(J[[n-1]]) + J[[n]]%*%(TPsr[[n]] - Phi%*%tPtt[[n]])%*%t(J[[n-1]]) } # Q(,) calculations st1 <- vector("list",N) St1 <- matrix(0,nrow=dim(TPsr[[N]])[1],ncol=dim(TPsr[[N]])[2]) rcov <- matrix(0,nrow=dim(R)[1],ncol=dim(R)[2]) for(n in seq(from=N,to=1)){ st1[[n]] <- TPsr[[n]] + Txs[[n+1]]%*%t(Txs[[n]]) St1 <- St1 + st1[[n]] rcov <- rcov + (y[,n]-A[[n]]%*%Txs[[n+1]])%*%t(y[,n]-A[[n]]%*%Txs[[n+1]]) + A[[n]]%*%TPss[[n+1]]%*%t(A[[n]]) } # update mu, Phi, Q, R # NB: structured parameter matrices can be used by fixing the appropriate elements mu <- Txs[[1]] Phi <- St1%*%solve(Ss0) Q <- matrix(1/N,nrow=dim(Q)[1],ncol=dim(Q)[2])*(St0 - St1%*%solve(Ss0)%*%t(St1)) R <- matrix(1/N,nrow=dim(R)[1],ncol=dim(R)[2])*rcov iter <- iter + 1 if(iter>1) cvg <- (-2*old.loglik) - (-2*loglik) old.loglik <- loglik } # end while loop r <- dim(R0)[1] BIC <- loglik - 0.5*(p + p^2 + r*(r+1)/2 + r*(r+1)/2)*log(N) # mu, Phi, Q, R parameters out <- list(mu,Phi,Q,R,BIC) names(out) <- c("mu","Phi","Q","R","BIC") out } forecast <- function(ahead,numobs,y,A,Sigma0,mu0,Phi0,Q0,R0){ # ahead - the desired number of forecasts of y[,t], t>T # y[,t] rx1 # A[[t]] rxp # Sigma0 pxp # mu0 px1 # Phi0 pxp # Q0 pxp # R0 rxr # EM estimation of Kalman filter # Reference: "An approach to time series smoothing # and forecasting using the EM algorithm." # R. H. Shumway and D. S. Stoffer # J. Time Series Anal. Vol. 3, No. 4 253-264 # y(t) = A(t)%*%x(t) + v(t), for t=1,...,T # x(t) = Phi%*%x(t-1) + w(t), for t=1,...,T # y rx1 observation vector # A rxp design matrix # E[v(t)]=0, R = E[v(t)%*%t(v(t))] rxr observation noise covariance # x px1 state vector # Phi pxp transition matrix # E[w(t)]=0, Q = E[w(t)%*%t(w(t))] pxp state noise covariance # initial vector x0, E[x0] = mu, Sigma = E[(x0-mu)%*%t(x0-mu)] # Notation: # sxt = E[x(t)|y(1),...,y(s)] # sPtu = E[(x(t)-sxt)%*%t(x(u)-sxu)|y(1),...,y(s)] N <- numobs + ahead mu <- mu0 Phi <- Phi0 Q <- Q0 R <- R0 r <- dim(A[[1]])[1] p <- dim(A[[1]])[2] # gain K pxr and covariances sPtt, tPtt sPtt <- vector("list",N) K <- vector("list",N) tPtt <- vector("list",N+1) tPtt[[1]] <- Sigma0 sxt <- vector("list",N) txt <- vector("list",N+1) txt[[1]] <- mu # forward recursions for(n in seq(N)){ sPtt[[n]] <- Phi%*%tPtt[[n]]%*%t(Phi) + Q # K(t)=0, for t>T if(n<=numobs) { K[[n]] <- sPtt[[n]]%*%t(A[[n]])%*%solve(A[[n]]%*%sPtt[[n]]%*%t(A[[n]])+R) } else { K[[n]] <- matrix(0,nrow=p,ncol=r) } tPtt[[n+1]] <- sPtt[[n]] - K[[n]]%*%A[[n]]%*%sPtt[[n]] # filter estimators sxt, txt sxt[[n]] <- Phi%*%txt[[n]] if(n<=numobs) { txt[[n+1]] <- sxt[[n]] + (K[[n]]%*%(y[,n]-(A[[n]]%*%sxt[[n]]))) } else { txt[[n+1]] <- sxt[[n]] } } # backward recursions J <- vector("list",N) Txs <- vector("list",N+1) Txs[[N+1]] <- txt[[N+1]] TPss <- vector("list",N+1) TPss[[N+1]] <- tPtt[[N+1]] # smoother estimators Txs for(n in seq(from=N,to=1)){ J[[n]] <- tPtt[[n]]%*%t(Phi)%*%solve(sPtt[[n]]) Txs[[n]] <- txt[[n]] + J[[n]]%*%(Txs[[n+1]]-sxt[[n]]) # covariance TPss TPss[[n]] <- tPtt[[n]] + J[[n]]%*%(TPss[[n+1]]-sPtt[[n]])%*%t(J[[n]]) } # covariance TPsr TPsr <- vector("list",N) TPsr[[N]] <- (diag(p) - K[[N]]%*%A[[N]])%*%Phi%*%tPtt[[N]] for(n in seq(from=N,to=2)){ TPsr[[n-1]] <- tPtt[[n]]%*%t(J[[n-1]]) + J[[n]]%*%(TPsr[[n]] - Phi%*%tPtt[[n]])%*%t(J[[n-1]]) } Ty <- matrix(0,nrow=r,ncol=N) Tx <- matrix(0,nrow=p,ncol=N) for(i in seq(from=N+1,to=2)){ Ty[,i-1] <- A[[i-1]]%*%Txs[[i]] Tx[,i-1] <- Txs[[i]] } out <- list(Ty,Tx) names(out) <- c("obs.pred","state.pred") out } (3) GENERAL DISCLAIMER This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 1, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. For a copy of the GNU General Public License write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. In short: you may use this code any way you like, as long as you don't charge money for it, remove this notice, or hold anyone liable for its results. Gary Nan Tie (2i2) 902-7265 nantig@fi.gs.com