#########################################
#  第６章 ARMAモデルによる時系列の解析
#########################################
#
# 白色雑音
w <- as.ts(rnorm(800))
plot(w,ylim=c(-3,3))
plot(as.ts(w[1:200]),ylim=c(-3,3))


# インパルス応答関数，自己共分散関数，PARCOR
# パワースペクトル，固有根
 a <- c(0.9*sqrt(3), -0.81)
 b <- c(0.9*sqrt(2), -0.81)

# AR model : y(n) = 0.9sqrt(3)y(n-1) - 0.81y(n-2) + v(n)
z <- armaimp(arcoef=a, v=1.0, n=1000, lag=20)
z$croot.ar

# MA model : y(n) = v(n) -0.9sqrt(2)v(n-1) + 0.81v(n-2)
z <- armaimp(macoef=b, v=1.0, n=1000, lag=20)
z$croot.ma

# ARMA model :  y(n) = 0.9sqrt(3)y(n-1) - 0.81y(n-2)
#                      + v(n) -0.9sqrt(2)v(n-1) + 0.81v(n-2)
z <- armaimp(arcoef=a, macoef=b, v=1.0, n=1000, lag=20)
z$croot.ar
z$croot.ma

# ラインスペクトル
# ARMA model :  y(n) = 0.99sqrt(2)y(n-1) - 0.99^2y(n-2)
#                      + v(n) -0.95sqrt(2)v(n-1) + 0.95^2v(n-2)
a <- c(0.99*sqrt(2), -0.9801)
b <- c(0.95*sqrt(2), -0.9025)
z <- armaimp(arcoef=a, macoef=b, v=1.0, n=1000, lag=20)
z$croot.ar
z$croot.ma

#####################
# 多変量ARモデル
#####################
# 
data(HAKUSAN) # Yaw rate, rolling, pitching and rudder angle for the ship on the open sea
length <- dim(HAKUSAN)[1] /2
y <- matrix(, length, 3) 
for( i in 1:length) { 
　　y[i,1] <- HAKUSAN[i*2,1]  # yaw rate
　　y[i,2] <- HAKUSAN[i*2,2]  # rolling 
　　y[i,3] <- HAKUSAN[i*2,4]} # rudder angle  
marfit(y,lag=20)
marlsq(y,lag=20)
　　z <- marfit(y,lag=20)