#########################################
#  第７章 ARモデルの推定
#########################################
# テストデータの読み込み
hakusan <- as.ts(read.csv("hakusan_new.csv"))
hakusan1 <- hakusan[,1]
sunspot <- as.ts(read.csv("sunspot_new.csv"))
maxtemp <- as.ts(read.csv("maxtemp.csv"))
blsfood <- as.ts(read.csv("blsfood_new.csv"))
whard <- as.ts(read.csv("whard_new.csv"))
mye1f <- as.ts(read.csv("mye1f_new.csv"))

# AR model fitting for Candian Lynx data
# method=1 (default) Yule-Walker method
# method=2 Householder least squares method
# method=3 Parcor method (Partial autoregrssion)
# method=4 PARCOR
# method=5 Burg's algorithm (MEM)
#
arfit(lynx)
arfit(lynx,method=2)
arfit(lynx,method=4)

# Other test data
arfit(hakusan1)
arfit(sunspot)
arfit(maxtemp)
arfit(blsfood)
arfit(whard)
arfit(mye1f)
arfit(mye1f,method=3)
arfit(mye1f,lag=20)



# MAR model fitting by Yule-Walker method for Hakusan data
data(HAKUSAN)   # Yaw rate, rolling, pitching and rudder angle for the ship on the open sea
  y <- as.matrix(HAKUSAN)
  z <- marfit(y,20)
  z$maice.order
  z$aic

plot(z$aic,type="h",ylim=c(46200,46500),lwd=8,col="blue")


# MAR model fitting by the Householder method
#
  data(HAKUSAN)   # Yaw rate, Rolling, pitching and rudder angle for the ship on the open sea
  length <- dim(HAKUSAN)[1]
  y <- matrix(, length, 3)
  y[,1] <- HAKUSAN[,1]
  y[,2] <- HAKUSAN[,2]
  y[,3] <- HAKUSAN[,4]
  z <- marlsq(y)
#
# Compute cross-spectra, coherency and relative noise contribution
#
  marspc(z$arcoef, z$v)