### 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)
}