#########################################
#  第9章 状態空間モデルによる欠測値の補間
#########################################
だた（BLSALLFOOD)
## AR model (l=1, k=1) AR次数は自動決定
#
## Example of interpolation of missing values (AR
model : m=15, k=1)
z2 <‐ arfit(BLSALLFOOD, plot = FALSE)
tau2 <‐ z2$sigma2
# m = maice.order, k=1
m2 <‐ z2$maice.order
arcoef <‐ z2$arcoef[[m2]]
f <‐ matrix(0.0e0, m2, m2)
f[1, ] <‐ arcoef
if (m2 != 1)
for (i in 2:m2) f[i, i‐1] <‐ 1
g <‐ c(1, rep(0.0e0, m2‐1))
h <‐ c(1, rep(0.0e0, m2‐1))
q <‐ tau2[m2+1]
r <‐ 0.0e0
x0 <‐ rep(0.0e0, m2)
v0 <‐ NULL
tsmooth(BLSALLFOOD, f, g, h, q, r, x0, v0, missed
= c(41, 101), np = c(30, 20))

###################################
# AR(5)を使う場合
#
# AR order=5
m2 <- ５
arcoef <- z2$arcoef[[m2]]
f <- matrix(0.0e0, m2, m2)
f[1, ] <- arcoef
if (m2 != 1)
for (i in 2:m2) f[i, i-1] <- 1
g <- c(1, rep(0.0e0, m2-1))
h <- c(1, rep(0.0e0, m2-1))
q <- tau2[m2+1]
r <- 0.0e0
x0 <- rep(0.0e0, m2)
v0 <- NULL
tsmooth(BLSALLFOOD, f, g, h, q, r, x0, v0,
missed = c(41, 101), np = c(30, 20))





#########################################
#  第10章 状態空間モデルによる時系列の解析
#########################################
#
# ARMA モデルの推定
#
data(Sunspot)         # Sun spot number data
y <- log10(Sunspot)   # y = log( sunspot data )

#
# 単にAR(4,2)を計算する場合
armafit(y, ar.order=4, ma.order=2)


#
# ARMA(6,3) modelを推定してスペクトル等を計算する場合
#
z1 <- armafit(y$z, ar.order=6, ma.order=3)
z1
# Estiamte ARMA(2,4) and compute Impuse response, autocovariance,
# PARCOR, power spectrum and characteristic roots.
len <- length(Sunspot)
z <- armaimp(arcoef=z1$arcoef,
macoef=z1$macoef, v=z1$sigma2, n=len,
lag=20)
armafit(y, ar.order=4, ma.order=2)


#
# ARMA(3,2) model　　#パラメータの初期値を指定する場合
# Use the parameters of ARMA(3,1) as the initial estimates
# of ARMA(3,2)
#
z1 <- armafit(y, ar.order=3, ma.order=2, ar=c(1.42758,-0.69843,0.00998), ma=c(0.3974,0.0))
z1



# 次数の全探索
# Automatic fitting of ARMA models
#
mmax <- 10 # maximum AR order
lmax <- 5 # maximum MA order
m2 <- (mmax+1)*(lmax+1)
aic <- matrix(1:m2,nrow=mmax+1,ncol=lmax+1)
aic[1,1] <- 0
ii <- seq(1,mmax,length=mmax)
jj <- seq(1,lmax,length=lmax)
#
for (i in ii){
z1 <- armafit(y$z,ar.order=i,ma.order=0) # AR(i) model
aic[i+1,1] <- z1$aic
for (j in jj){
ar <- z1$arcoef
ma <- z1$macoef
ma[j] <- 0.001
z1 <- armafit(y$z,ar.order=i,ma.order=j,ar,ma) # ARMA(i,j)
aic[i+1,j+1] <- z1$aic
}
}
for (j in jj){ # MA(j) models (j=1,…,lmax)
ma <- z1$macoef
z1 <- armafit(y$z,ar.order=0,ma.order=j,ma)
aic[1,j+1] <- z1$aic
}


####################################
# 第11章　トレンドの推定
####################################
#
# Polinomial regression model 
#   highest order = 7
#
  data(Temperature)   # Highest Temperature Data of Tokyo
  polreg(Temperature, 7)
#
# plot AIC's of polinomial regression models
#
plot(aic,type="b",lwd=2,ylim=c(2450,2520))

#
#  Trigonometric regression models
# 
 lsqr(Temperature)
#

####################################
#  Trend Model
####################################
#
#  Second order trend model
#
trend(Temperature, trend.order =2)

#
#  First order trend model
#     TAU2 = 0.223, Sig2 = 0.001
#
 trend(Temperature, trend.order = 1, tau2.ini = 0.223, delta = 0.001)