### Selected R commands and output used in Chapter 6 require(dlm) #### Build function to MLE estimation buildTrackLinear=function(parm){ q <- parm[1:3] r <- parm[4:6] H <- c(1,0,0,1,0,0, 0,1,0,0,1,0, 0,0,1,0,0,1, 0,0,0,1,0,0, 0,0,0,0,1,0, 0,0,0,0,0,1) H <- matrix(H,ncol=6,nrow=6,byrow=T) G <- c(1,0,0,0,0,0, 0,1,0,0,0,0, 0,0,1,0,0,0) G <- matrix(G,ncol=6,nrow=3,byrow=T) W <- c(0.5*q[1], 0, 0, 0, 0.5*q[2], 0, 0, 0, 0.5*q[3], q[1], 0, 0, 0, q[2], 0, 0, 0, q[3]) W <- matrix(W,ncol=3,nrow=6,byrow=T) V <- diag(r) return(list( m0=rep(0,6), C0=1000000*diag(6), FF=G, GG=H, V=V%*%t(V), W=W%*%t(W) )) } #### end function ## #true parameters q=c(0.1,0.1,0.02) r <- c(1,1,1) ### Run Kalman filter with DLM package library(dlm) myModel <- buildTrackLinear(c(q,r)) #true model outFilter <- dlmFilter(y,myModel) outSmooth <- dlmSmooth(y,myModel) ###Estimation par.init <- c(1,1,1,10,10,10) outMLE <- dlmMLE(y,par.init,buildTrackLinear) #### library(dlm) #### Build function to MLE estimation: buildVaryingBeta1=function(parm){ phiW <- parm[1:5] phiV <- parm[6] H <- diag(c(1,1,1,1,1)) W0 <- diag(phiW) W0 <- W0%*%t(W0) W <- matrix(W0,ncol=5,nrow=5,byrow=T) G <- c(1,1,1,1,1) G <- matrix(G,ncol=5,nrow=1,byrow=T) JGG <- c(0,1:4) #JGG indicates which element is time varying JGG <- matrix(JGG,ncol=5,nrow=1,byrow=T) V <- matrix(phiV[1]**2,ncol=1,nrow=1) TXX <- xx[,1:4] #TXX provides the time varying inputs return(list(m0=mm0,C0=CC0,FF=G,GG=H,V=V,W=W,JFF=JGG,X=TXX)) } ### end function ### xxx <- read.csv(file="FF_example.csv") y <- xxx[,2] xx <- xxx[,3:6] ## FF factor plot year0 <- 1946:2010 y <- xxx[,2] xx <- xxx[,3:6] name <- c('cnst','market','SMB','HML','MoM') par(mfrow=c(2,2),mai=0.7*c(1,1,1,1)) for(j in 1:4){ plot(year0,xx[(year0-1945)*12,j],type='l',xlab='time',ylab='return',main=name[j-1]) abline(h=0,lty=2) } ## fig_varying_beta0.pdf ## FF factor plot ## y-plot par(mfrow=c(1,1),mai=1*c(1,1,1,1)) plot(year0,x[(year0-1945)*12,1],type='l',xlab='time',ylab='return',main='portfolio') # overall FF model out0 <- lm(y~xx) summary(out0) #rolling window estimation ntotal <- length(y) nyear <- ntotal/12 beta <- matrix(ncol=5,nrow=(nyear-2)) for(i in 3:nyear){ t0 <- (i-3)*12+(1:36) y0 <- y[t0] ##excess return of portfolio 22 xx0 <- x[t0,] out <- lm(y0~xx0) beta[i-2,] <- out\$coef } year <- 1948:2010 par(mfrow=c(2,2),mai=0.7*c(1,1,1,1)) for(j in 2:5){ ylim <- range(cbind(ucl[,j],lcl[,j])) plot(year,beta[,j],type='l',xlab='time',ylab='beta',main=name[j-1],ylim=ylim) lines(year,year*0+summary(out0)\$coef[j,1]) lines(year,year*0+summary(out0)\$coef[j,1]+summary(out0)\$coef[j,2]*2,lty=2) lines(year,year*0+summary(out0)\$coef[j,1]-summary(out0)\$coef[j,2]*2,lty=2) lines(year,ucl[,j],lty=2) lines(year,lcl[,j],lty=2) } ## fig_varying_beta1.pdf ## 3-year regression ##++++++ use diffused initial value mm0 <- c(1,1,1,1,1) CC0 <- diag(c(1,1,1,1,1))*1000 ww <- c(1,1,1,1,1) vv <- 1 init.par <- c(ww,vv) MLE2 <- dlmMLE(y,init.par,buildVaryingBeta1,debug=F) MLE2 MyModel2 <- buildVaryingBeta1(MLE2\$par) ss2 <- dlmSmooth(y,MyModel2) #--- smoothed state plot year <- 1948:2010 par(mfrow=c(2,2),mai=0.7*c(1,1,1,1)) for (j in 2:5){ plot(year,beta[,j],type='p',xlab='time',ylab='beta',main=name[j-1]) lines(year,ss2\$s[(3:65)*12+1,j]) } ##++++++ use an arbitrary state initial value mm0 <- c(1,2,-0.2,0.8,-0.3) CC0 <- diag(c(0.1,0.01,0.1,10,10)) ww <- c(30,30,30,10,10) vv <- 1 init.par <- c(ww,vv) MLE <- dlmMLE(y,init.par,buildVaryingBeta1,debug=F) MyModel <- buildVaryingBeta1(MLE\$par) ss <- dlmSmooth(y,MyModel) # year <- 1948:2010 par(mfrow=c(2,2),mai=0.7*c(1,1,1,1)) for (j in 2:5){ plot(year,beta[,j],type='p',xlab='time',ylab='beta',main=name[j-1]) lines(year,ss\$s[(3:65)*12+1,j]) } #### ### Build function to MLE estimation: sequential processing buildMixedFreqSeq4=function(parm){ nyear <- 30 phiH <- parm[1:10] phiW <- parm[11:20] phiV <- parm[21:24] H <- c(phiH[1],0,0,phiH[2],0,0, 0,phiH[3],0,phiH[4],0,0, 0,0,phiH[5],phiH[6],0,0, phiH[7:10],0,0, 0,0,0,1,0,0, 0,0,0,0,1,0) H <- matrix(H,ncol=6,nrow=6,byrow=T) JHH <- c(1,0,0,2,0,0, 0,3,0,4,0,0, 0,0,5,6,0,0, 7:10,0,0, 0,0,0,11,12,0, 0,0,0,0,13,14) JHH <- matrix(JHH,ncol=6,nrow=6,byrow=T) W0 <- matrix(c(phiW[1:4], 0,phiW[5:7], 0,0,phiW[8:9], 0,0,0,phiW[10]),ncol=4,nrow=4,byrow=T) W0 <- W0%*%t(W0) W1 <- rep(0,16) W <- c(c(W0[1,]),0,0, c(W0[2,]),0,0, c(W0[3,]),0,0, c(W0[4,]),0,0, rep(0,6), rep(0,6)) W <- matrix(W,ncol=6,nrow=6,byrow=T) JWW <- c(15:18,0,0, 19:22,0,0, 23:26,0,0, 27:30,0,0, rep(0,6), rep(0,6)) JWW <- matrix(JWW,ncol=6,nrow=6,byrow=T) G <- c(1,0,0,0,0,0) G <- matrix(G,ncol=6,nrow=1,byrow=T) JGG <- c(31:36) JGG <- matrix(JGG,ncol=6,nrow=1,byrow=T) V <- matrix(phiV[1]**2,ncol=1,nrow=1) JVV <- matrix(37,ncol=1,nrow=1) xx1 <- c(phiH[1:10],1,0,1,0, c(W0), 1,rep(0,5), phiV[1]**2) xx2 <- c(1,0,1,0,1,0,0,0,0,1,0,1,0,1, c(W1), 0,1,rep(0,4), phiV[2]**2) xx3 <- c(1,0,1,0,1,0,0,0,0,1,0,1,0,1, c(W1), 0,0,1,rep(0,3), phiV[3]**2) xx4 <- c(1,0,1,0,1,0,0,0,0,1,0,1,0,1, c(W1), 0,0,0,1,1,1,phiV[4]**2) xx <- c(rep(c(xx1,xx2,xx3),3),xx4) TXX <- matrix(rep(xx,nyear),ncol=37,nrow=nyear*40,byrow=T) return(list( m0=rep(0,6), C0=1000000*diag(6), FF=G, GG=H, V=V, W=W, JFF=JGG, JGG=JHH, JWW=JWW, JVV=JVV, X=TXX )) } ### end function ### # Analysis y_start=22 ## 1948+22=1970 nyear <- 30 # temp <- c(t(cbind(CPIuse_m,INDPuse_m,UEMPuse_m))) yseq <- 1:(nyear*40) for(i in 1:(nyear*4)){ yseq[(i-1)*10+(1:9)] <- temp[(i-1)*9+(1:9)] yseq[(i-1)*10+10] <- GDPuse_m[i] } # phi <- c(0.7,-0.3,0.7,-0.3,0.7,-0.3,0.1,0.1,0.1,0.8) ww <- c(0.3,0,0,0,0.5,0,0,0.5,0,0.3) vv <- c(1,1,1,1) parm <- c(phi,ww,vv) MLE <- dlmMLE(yseq,parm,buildMixedFreqSeq4) # MyModel <- buildMixedFreqSeq4(MLE\$par) ss <- dlmSmooth(yseq,MyModel) ## ss provides smoothed states #### ####### buildDFM3 buildDFM3=function(parm){ hh <- c(parm[1],parm[2],0,parm[3]) gg <- c(1,1,parm[4:11]) qq <- c(parm[12:13],0,parm[14]) rr <- parm[15:19] H <- matrix(hh,ncol=2,nrow=2,byrow=T) G <- matrix(gg,ncol=2,nrow=5,byrow=T) W <- matrix(qq,ncol=2,nrow=2,byrow=T) V <- diag(rr) return(list( m0=rep(0,2), C0=5*diag(2), FF=G, GG=H, V=V%*%t(V), W=W%*%t(W) )) } ### end function ### ### DFM outsample prediction outsample.prediction.DFM=function(x,np,n.Start,n.End,par.init){ xpred <- x SS <- rep(0,5) for(i in (n.Start-np):(n.End-np)){ MLE <- dlmMLE(x[1:i,],par.init,buildDFM3) myMode <- buildDFM3(MLE\$par) filter <- dlmFilter(x[1:i,],myMode) pp <- dlmForecast(filter,nAhead=np) xpred[i+np,] <- pp\$f[np,] SS <- SS+(x[i+np,]-xpred[i+np,])**2 } return(list(xpred=xpred,SS=SS,meanSS=mean(SS))) } ### PCA outsample prediction outsample.prediction.PCA=function(x,np,n.Start,n.End){ xpred <- x SS <- rep(0,5) for(i in (n.Start-np):(n.End-np)){ PCA <- princomp(x[1:i,], cor=F) ff <- PCA\$score[,1:2] out <- VAR(ff) pp <- predict(out,n.ahead=np) pp0 <- c(pp\$fcst\$Comp.1[np,1],pp\$fcst\$Comp.2[np,1]) xpred[i+np,] <- PCA\$loading[,1:2]%*%pp0 SS <- SS+(x[i+np,]-xpred[i+np,])**2 } return(list(xpred=xpred,SS=SS,meanSS=mean(SS))) } ### univariate AR prediction outsample.prediction.uAR=function(x,np,n.Start,n.End){ xpred <- x SS <- rep(0,5) for(i in (n.Start-np):(n.End-np)){ for(j in 1:5){ out <- arima(x[1:i,j],c(1,0,0)) xpred[i+np,j] <- predict(out,n.ahead=np)\$pred[np] } SS <- SS+(x[i+np,]-xpred[i+np,])**2 } return(list(xpred=xpred,SS=SS,meanSS=mean(SS))) } ### Analysis ### xx.nm <- matrix(scan("GDP_panel_nm.dat"),ncol=5) #-------- PCA analysis ---------------# PCAout <- princomp(xx.nm, cor=F) ff <- PCAout\$scores[,1:2] library(vars) out <- VAR(ff) #-------- DFM analysis ---------------# hh3 <- c(0.7,0.2,0.8) gg3 <- rep(c(0.7,0.7), 4) qq3 <- c(1,-1,1) rr3 <- rep(0.3,5) parm3 <- c(hh3,gg3,qq3,rr3) # initial parameter values outMLE3 <- dlmMLE(xx.nm,parm3,buildDFM3) myModel3 <- buildDFM3(outMLE3\$par) outSmooth3 <- dlmSmooth(xx.nm,myModel3) fit3 <- xx.nm for(i in 1:5){ fit3[,i] <- outSmooth3\$s[-1,]%*%myModel3\$F[i,] } print(sum(xx.nm-fit3)**2)/sum(xx.nm**2)) #---------- Prediction ---------------# for(np in 1:3){ pred.d <- outsample.prediction.DFM(xx.nm,np,88,107,outMLE3\$par) pred.p <- outsample.prediction.PCA(xx.nm,np,88,107) pred.ar <- outsample.prediction.uAR(xx.nm,np,88,107) }