SARIMAとかで遊んでみたメモ

完全に自分用。

x <- 1:20
y <- c()
for(i in x){
 y[i] <- 3*i^2+2
}
y

y <- as.ts(y)
#一回の差分を取る
diff(y)
#n回の差分を取るためにはdifferencesという引数を取る
diff(y,differences=3)

ibm <- c(460,457,452,459,462,459,463,479,493,490,492,498,499,497,496,490,489,487,487,491,487,482,479)
ibm

ts.plot(diff(ibm,differences=3))

arima(ibm,order=c(2,3,1))
plot(ibm)
acf(ibm)

library(MASS)
acf(deaths)
plot(deaths)

mean(diff(deaths,differences=12))
?arima.sim
example(arima.sim)

ts.sim <- arima.sim(list(order=c(1,1,0),ar=0.7),n=200)
plot(ts.sim)
plot(diff(ts.sim))
acf(diff(ts.sim))

acf(ts.sim)


ts.sim <- arima.sim(list(order=c(1,0,0),ar=0.9),n=200)
plot(ts.sim)
lines(lowess(ts.sim))
plot(diff(ts.sim))
acf(diff(ts.sim))

diff(deaths,lag=12,differences=3)

deaths.diff <- diff(deaths,12)
deaths.diff
acf(deaths,30)
acf(deaths.diff,30)
acf(deaths.diff,30,type="partial")


arima(deaths,order=c(1,0,1))
arima(deaths,order=c(1,1,1))


co2
plot(co2)
#季節性を考慮することで線形関係のトレンドを取り除く
diff(co2,lag=12)
plot(diff(co2,lag=12))
#二回目の差分を取ることでレベルの違いを取り除く
diff(co2,lag=12,differences=2)
plot(diff(co2,lag=12,differences=2))

diff(co2,differences=2)
plot(diff(co2,differences=2))
acf(diff(co2,differences=2))


data(sunspot)
sunspot.month
plot(sunspot.month, xlab="", ylab="sunspot")
acf(sunspot.month,  xlab="", main="")
plot(diff(sunspot.month),xlab="", ylab="diff(sunspot)")
acf(diff(sunspot.month), xlab="", main="")

plot(diff(sunspot.month,lag=24))


my.sarima.sim <- function (
    n = 20, 
    period = 12, 
    model, 
    seasonal
) {    
  x <- arima.sim( model, n*period )
  x <- x[1:(n*period)]
  for (i in 1:period) {
    xx <- arima.sim( seasonal, n )
    xx <- xx[1:n]
    x[i + period * 0:(n-1)] <- 
      x[i + period * 0:(n-1)] + xx
  }
  x <- ts(x, frequency=period)
  x
}

x <- my.sarima.sim(
  20, 
  12, 
  list(ar=.6, ma=.3, order=c(1,0,1)),
  list(ar=c(.5), ma=c(1,2), order=c(1,1,2))
)

x
plot(x)
#ちっとも低減する兆しがない
acf(x)
#周期性を考慮できていない
acf(diff(x))
#sinカーブ減衰っぽくなった
#季節のラグを考えてから、差分を一回取った
acf(diff(x,lag=12))

x <- my.sarima.sim(
  20, 
  12, 
  list(ar=.6, ma=.3, order=c(1,1,1)),
  list(ar=c(.5), ma=c(1,2), order=c(1,1,2))
)

plot(x)
#周期性の影響が強く出ている。低減しない。
acf(x,50)
#非定常なデータになっている
acf(diff(x,lag=12))
#こっちだと季節性で差分を取ってしまっている
acf(diff(x,lag=12,differences=2))

acf(diff(diff(x,lag=12)))