#########################################
#  第9章 状態空間モデルによる時系列の解析
#########################################
#
# AR モデルによる長期予測
#
## AR model (l=1, k=1)  長期予測
#
  data(Blsallfood)
  z1 <- arfit(Blsallfood, plot=FALSE)
  m <- z1$maice.order
  tau2 <- z1$sigma2[m+1]
  arcoef <- z1$arcoef

  f <- matrix(0.0e0, m, m)
  f[1,] <- arcoef
  for( i in 2:m ) f[i,i-1] <- 1
  g <- c(1, rep(0.0e0, m-1))
  h <- c(1, rep(0.0e0, m-1))
  q <- tau2
  r <- 0.0e0
  x0 <- rep(0.0e0, m)
  v0 <- NULL

  z <- tsmooth(Blsallfood, f, g, h, q, r, x0, v0, 156, 156,
               missed=c(120), np=c(36))
# plot mean vector and estimation error
  xss <- z$mean.smooth[1,] + mean(Blsallfood)
  cov <- z$cov.smooth
  c1 <- xss + sqrt(cov[1,])
  c2 <- xss - sqrt(cov[1,])
  err <- z$esterr
  par(mfcol=c(2,1))
  ymax <- as.integer(max(xss,c1,c2)+1)
  ymin <- as.integer(min(xss,c1,c2)-1)
  plot(c1, type='l', ylim=c(ymin,ymax), col=2,
       xlab="Mean vectors of the smoother XSS(1,) +/- standard deviation",
       ylab="")

  par(new=TRUE)
  plot(c2, type='l', ylim=c(ymin,ymax), col=3, xlab="", ylab="")

  par(new=TRUE)
  plot(xss, type='l', ylim=c(ymin,ymax), xlab="", ylab="")
  plot(err[,1,1], type='h', xlim=c(1,length(xss)), xlab="estimation error",
       ylab="")



######################################
#
#  AR次数 1　（下記のものに置き換える）
#
z1 <- arfit(Blsallfood, plot=FALSE,lag=1)
#
#  AR次数 5　（下記のものに置き換える）
#
z1 <- arfit(Blsallfood, plot=FALSE,lag=5)


###############################################################
#
# ARモデルによる欠測値の補間
#
## AR model (l=1, k=1)　　AR次数は自動決定
#
data(Blsallfood) z1 <- arfit(Blsallfood, plot=FALSE) 
m <- z1$maice.order 
tau2 <- z1$sigma2[m+1] 
arcoef <- z1$arcoef 

f <- matrix(0.0e0, m, m) 
f[1,] <- arcoef 
for( i in 2:m ) f[i,i-1] <- 1 
g <- c(1, rep(0.0e0, m-1)) 
h <- c(1, rep(0.0e0, m-1)) 
q <- tau2 r <- 0.0e0 
x0 <- rep(0.0e0, m) 
v0 <- NULL 
z <- tsmooth(Blsallfood, f, g, h, q, r, x0, v0, 156, 170, missed=c(41,101), np=c(30,20)) 
# plot mean vector and estimation error 
xss <- z$mean.smooth[1,] + mean(Blsallfood) 
cov <- z$cov.smooth 
c1 <- xss + sqrt(cov[1,]) 
c2 <- xss - sqrt(cov[1,]) 
err <- z$esterr 
par(mfcol=c(2,1)) 
ymax <- as.integer(max(xss,c1,c2)+1) 
ymin <- as.integer(min(xss,c1,c2)-1) 
plot(c1, type=‘l’, ylim=c(ymin,ymax), col=2, xlab=“Mean 　　vectors of the smoother XSS(1,) +/- standard deviation", ylab="") 
par(new=TRUE) 
plot(c2, type='l', ylim=c(ymin,ymax), col=3, xlab="", ylab="") 
par(new=TRUE) 
plot(xss, type='l', ylim=c(ymin,ymax), xlab="", ylab="") plot(err[,1,1], type='h', xlim=c(1,length(xss)), xlab="estimation error", ylab="") 


# plot mean vector and estimation error
  xss <- z$mean.smooth[1,] + mean(Blsallfood)
  cov <- z$cov.smooth
  c1 <- xss + sqrt(cov[1,])
  c2 <- xss - sqrt(cov[1,])
  err <- z$esterr
  par(mfcol=c(2,1))
  ymax <- as.integer(max(xss,c1,c2)+1)
  ymin <- as.integer(min(xss,c1,c2)-1)
  plot(c1, type='l', ylim=c(ymin,ymax), col=2,
       xlab="Mean vectors of the smoother XSS(1,) +/- standard deviation",
       ylab="")

  par(new=TRUE)
  plot(c2, type='l', ylim=c(ymin,ymax), col=3, xlab="", ylab="")

  par(new=TRUE)
  plot(xss, type='l', ylim=c(ymin,ymax), xlab="", ylab="")
  plot(err[,1,1], type='h', xlim=c(1,length(xss)), xlab="estimation error",
       ylab="")

#########################
#
#  AR次数 1　（下記のものに置き換える）
#
z1 <- arfit(Blsallfood, plot=FALSE,lag=1)
#
#  AR次数 5
#
z1 <- arfit(Blsallfood, plot=FALSE,lag=5)


######################################
## Trend model (l=3, m=2, k=2)
#  n <- 400  or
  n <- 500
  l <- 3
  m <- 2
  k <- 2
  f <- matrix(c(1, 0,  0, 1), m, m, byrow=TRUE)
  g <- matrix(c(1, 0,  0, 1), m, k, byrow=TRUE)
  h <- matrix(c(0.1, -0.1,  -0.05, 0.05,  0.2, 0.15), l, m, byrow=TRUE)
  q <- matrix(c(0.2*0.2, 0,  0, 0.3*0.3), k, k, byrow=TRUE)
  r <- matrix(c(0.2*0.2, 0, 0,  0, 0.1*0.1, 0,  0, 0, 0.15*0.15), l, l,
              byrow=TRUE)
  Xn <- matrix(0, m, n)
  x1 <- rnorm(n+100, 0, 0.2)
  x2 <- rnorm(n+100, 0, 0.3)
  x1 <- cumsum(x1)[101:(n+100)]
  x2 <- cumsum(x2)[101:(n+100)]
  Xn[1,] <- x1-mean(x1)
  Xn[2,] <- x2-mean(x2)
  Yn <- matrix(0, l, n)
  Wn <- matrix(0, l, n)
  Wn[1,] <- rnorm(n, 0, 0.2)
  Wn[2,] <- rnorm(n, 0, 0.1)
  Wn[3,] <- rnorm(n, 0, 0.15)
  Yn <- h %*% Xn + Wn
  Yn <- aperm(Yn, c(2,1))
  x0 <- c(Xn[1,1], Xn[2,1])
  v0 <- matrix(c(var(Yn[,1]), 0,  0, var(Yn[,2])), 2, 2, byrow=TRUE)
  npe <- n+20
  z <- tsmooth(Yn, f, g, h, q, r, x0, v0, n, npe, missed=n/2, np=n/20)
 
# plot mean vector and state vector
  xss <- z$mean.smooth
  par(mfcol=c(m,1))
  for( i in 1:m ) {
    ymax <- as.integer(max(xss[i,],Xn[i,])+1)
    ymin <- as.integer(min(xss[i,],Xn[i,])-1)
    plot(Xn[i,], type='l', xlim=c(1,npe), ylim=c(ymin,ymax),
         xlab=paste(" red : mean.smooth[",i,",]  /  black : Xn[",i,",]"),
         ylab="")
    par(new=TRUE)
    plot(xss[i,], type='l', ylim=c(ymin,ymax), xlab="", ylab="", col=2)
  }