Untitled

基礎知識

ARIMAモデルというのは、時系列値列がMoving average(ma)という要素とAutoregression integrated(ar)という要素とを想定した確率過程

• maはfilter()関数による“convolution”的フィルタリングが働いていると考えるモデル – これは、数値列が局所平均があって、それが数値列の位置によって変動するというモデル
• raはfilter()関数による“recursive”的フィルタリングが働いていると考えるモデル – これは、未来が現在と過去の値に上乗せする形で決まると考えるモデル

http://wolfweb.unr.edu/~zal/STAT758/Lab_3.pdf のAppendixにarima.sim()を使った実習資料があるので、それに沿って、ARIMAモデルに沿ったデータ作成をしながら、ARIMAモデルについて理解していくことにする

正規乱数とフィルタリングの基礎

基本は正規乱数

Z1<-rnorm(100)

変化させないフィルタリング、単純累積和を取るフィルタリング

フィルタをかける。convolutionフィルタとrecursiveフィルタの2つがある。 この２つのフィルタタイプは、ARIMA的にはautoregressive integratedとmoving averageとに相当するので、arフィルタ、maフィルタとも呼ぶことにする。

両者はフィルタとして何を用いるかで変化するが、最も基本的なフィルタを使うと、超基本的操作ができる。 maでは、何も変えない、という操作。

f1 <- c(1)
X1.ma <- filter(Z1,f1,method="convolution")
plot(Z1,X1.ma)

arでは、累積和、という操作。

f1 <- c(1)
X1.ar <- filter(Z1,f1,method="recursive")
plot(cumsum(Z1),X1.ar)

平滑化

maフィルタを掛けて平滑化する

次のようなフィルタを使う。指定数n=1,2,…に対して、長さ2n+1のベクトルを用意し、要素の値が、標準正規分布確率密度関数値(-n),(-n+1),…,0,1,…,nに比例した重みを割り付けるようなフィルタ

フィルタの重みベクトルは、自身とその前後に掛かる

N は2n+1より大きい数。 黒点が乱数。 赤線はmaフィルタ。

X <- rnorm(10000)
inputPanel(
  sliderInput("N", label = "no values:",
              min = 1, max = 5000, value = 1000, step = 1),
  sliderInput("n", label = "value specifies length of filter:",
              min = 0, max = 1000, value = 0, step = 1)
  )
  
renderPlot({
  n <- input$n
  #tmp <- dnorm((-n):n,0,input$r)
  tmp <- rep(1,2*n+1)
  fil <- tmp/sum(tmp)
  
  x <- X[1:input$N]
  fout <- filter(x,fil,method="convolution")
  #fout2 <- filter(x,fil,method="recursive")
  #ylim <- range(na.omit(c(x,cumsum(x),fout,fout2)))
  ylim <- range(na.omit(c(x,fout)))
  yr <- ylim[2]-ylim[1]
  ylim[1] <- ylim[1] - yr*0.1
  ylim[2] <- ylim[2] + yr*0.1
  plot(x,pch=20,ylim=ylim)
  points(fout,type="l",col=2,lwd=5)
  #points(cumsum(x),pch=1,col=1)
  #points(fout2,type="l",col=3)

  
})

再帰フィルタ(Recursive filter)

再帰フィルタは、過去の情報を累積的に、どんどんため込む形をしている。 したがって、単調に増加していくときには、指数関数的な増加が現れる。

逆に、単純な累積和であると、増えたり減ったりが激しく見えても、一貫して増えているわけではなく、増えたり減ったりの繰り返しであるなら、平坦な値列に変換するようなフィルタである。

フィルタベクトルは、自身より前の値に掛かる。

x <- rep(1,10)
y1 <- filter(x,c(1),method="recursive")
y2 <- filter(x,c(1,1),method="r")
y3 <- filter(x,c(1,1,1),method="r")
matplot(cbind(y1,y2,y3),type="l")

具体的にコードがあった方がわかるので、作ってみると、頭の方から、filter後の値を作り、それを累加していくようになっていることがわかる。

my.recursive.filter <- function(x,f){
  ret <- x
    for(i in 1:length(x)){
        for(j in 1:length(f)){
            if(i-j>0){
                ret[i] <- ret[i] + f[j]*ret[i-j]
            }
        }
    }
    ret
}

x <- rnorm(10)
my.recursive.filter(x,c(1,1))

##  [1] -0.2404171  1.0264624  0.6695158  2.6553075  3.6595878  7.1253772
##  [7] 11.3306302 18.6408791 31.2770960 49.5714631

filter(x,c(1,1),"r")

## Time Series:
## Start = 1 
## End = 10 
## Frequency = 1 
##  [1] -0.2404171  1.0264624  0.6695158  2.6553075  3.6595878  7.1253772
##  [7] 11.3306302 18.6408791 31.2770960 49.5714631

N はn+1より大きい数。 黒点が乱数の単純累積和。 赤線はraフィルタ。

X <- rnorm(10000)
inputPanel(
  sliderInput("N2", label = "no values:",
              min = 1, max = 5000, value = 1000, step = 1),
  sliderInput("n2", label = "value specifies length of filter:",
              min = 0, max = 10, value = 0, step = 1)
  )
  
renderPlot({
  n <- input$n2
  #tmp <- dnorm((-n):n,0,input$r)
  tmp <- rep(1,n+1)
  fil <- tmp/sum(tmp)
  
  x <- X[1:input$N2]
  #fout <- filter(x,fil,method="convolution")
  fout2 <- filter(x,fil,method="recursive")
  #ylim <- range(na.omit(c(x,cumsum(x),fout,fout2)))
  ylim <- range(na.omit(c(cumsum(x),fout2)))
  yr <- ylim[2]-ylim[1]
  ylim[1] <- ylim[1] - yr*0.1
  ylim[2] <- ylim[2] + yr*0.1
  plot(cumsum(x),pch=20,ylim=ylim)
  points(fout2,type="l",col=2,lwd=5)

})

フィルタを使って時系列データをシミュレーション作成する

正規乱数列を発生し、それをそのまま、値列とするとき、maフィルタを使うなら

x <- rnorm(10)
x.ma <- filter(x,c(1),method="convolution",sides = 1L) # sides=1Lは過去側のみのフィルタリングを指定
x

##  [1]  0.3253160 -1.0420168 -1.3197773  0.2885019 -0.7989580  1.2138630
##  [7] -1.1786703 -0.5989648  0.4553541 -0.1078745

x.ma

## Time Series:
## Start = 1 
## End = 10 
## Frequency = 1 
##  [1]  0.3253160 -1.0420168 -1.3197773  0.2885019 -0.7989580  1.2138630
##  [7] -1.1786703 -0.5989648  0.4553541 -0.1078745

plot(x.ma)

raモデルを使うなら

x.ra <- filter(x,c(0),method="recursive")
x

##  [1]  0.3253160 -1.0420168 -1.3197773  0.2885019 -0.7989580  1.2138630
##  [7] -1.1786703 -0.5989648  0.4553541 -0.1078745

x.ra

## Time Series:
## Start = 1 
## End = 10 
## Frequency = 1 
##  [1]  0.3253160 -1.0420168 -1.3197773  0.2885019 -0.7989580  1.2138630
##  [7] -1.1786703 -0.5989648  0.4553541 -0.1078745

plot(x.ra)

フィルタの長さを固定して、フィルタを色々に変えてどのようになるかを見てみよう。 ma = (1,0,0)の場合とar = (0,0,0)の場合は、発生した乱数そのもので、白色ノイズとなっている。

maの方は３パラメタを動かしてもがたぼこの具合に若干の違いが出るだけであるが、arの方は、急に値が発散したりする様子がわかる。

また、ar = (1,0,0)のときは、単純な酔歩になる様子も見て取れる。

X <- rnorm(10000)
inputPanel(

  sliderInput("N3", label = "no values:",
              min = 1, max = 5000, value = 1000, step = 1),
  sliderInput("ar1", label = "ar1:",
              min = -2, max = 2, value = 0, step = 0.01),
  sliderInput("ar2", label = "ar2:",
              min = -2, max = 2, value = 0, step = 0.01),
  sliderInput("ar3", label = "ar3:",
              min = -2, max = 2, value = 0, step = 0.01),
  sliderInput("ma1", label = "ma1:",
              min = -2, max = 2, value = 1, step = 0.01),
  sliderInput("ma2", label = "ma2:",
              min = -2, max = 2, value = 0, step = 0.01),
  sliderInput("ma3", label = "ma3:",
              min = -2, max = 2, value = 0, step = 0.01)
    
  )
  
renderPlot({
  N <- input$N3
  ar <- c(input$ar1,input$ar2,input$ar3)
  ma <- c(input$ma1,input$ma2,input$ma3)
  x <- X[1:N]
  foutma <- filter(x,ma,method="convolution",sides=1L) # sides=1Lは過去側のみのフィルタリングを指定
  foutar <- filter(x,ar,method="recursive")

  #ylim <- range(na.omit(c(x,cumsum(x),fout,fout2)))
  ylim <- range(na.omit(c(foutma,foutar)))
  yr <- ylim[2]-ylim[1]
  ylim[1] <- ylim[1] - yr*0.1
  ylim[2] <- ylim[2] + yr*0.1
  plot(foutma,type="l",ylim=ylim)
  points(foutar,type="l",col=2)

})

ARIMAモデル

RのARIMAモデル関数

ARIMAモデルはmaフィルタで説明できる部分とarフィルタで説明できる部分とが合わさって時系列データを作るモデルであり、 実データに対してARIMAモデルを適用し、そのmaフィルタ・arフィルタ等を推定する。

RのARIMA関係関数には、その推定をするarima()関数とARIMAモデルに基づいて時系列データをシミュレーション作成する関数arima.sim()がある。

コードは以下で表示できる。

# arima
# arima.sim

arima.sim() 関数の引数

まず、arima.sim()の引数を確認しよう

args(arima.sim)

## function (model, n, rand.gen = rnorm, innov = rand.gen(n, ...), 
##     n.start = NA, start.innov = rand.gen(n.start, ...), ...) 
## NULL

model 引数

引数 modelはもっとも基幹になる引数で

model <- list(ar = c(0.3,0.2),ma = c(0.1,0.1))

のように2ベクトルからなるリストである。 これが関数内部で以下のように使われる

if (length(model$ma)) {
        x <- filter(x, c(1, model$ma), sides = 1L) # method = "convolution" がデフォルト
        x[seq_along(model$ma)] <- 0 # NAを0にする処理
    }
    if (length(model$ar)) 
        x <- filter(x, model$ar, method = "recursive")

初めにmaフィルタを掛け、次にarフィルタを掛けている。 maフィルタのsides=1Lとは、フィルタを過去側のみに作用させる手法である。

また、maフィルタの冒頭に1が付加されている。これは、maがなくarもないときに、maフィルタがc(1)となるようにするための処理である。

したがって、arima.sim()関数の骨格は、ma型とar型のフィルタリングの連続実施であることがわかる。

model$order 引数

それ以外の処理は、“stop(”‘model’ must be list“)”のように、処理を中止する条件確認がいくつかあるほかは、次の処理が大きな処理である。

x <- diffinv(x, differences = d)

これは、ma,raフィルタリングの結果であるxを引数dで処理する行である。 このdは基幹引数であるmodelというリストのオプショナルな1成分で、長さ3のベクトルの第二要素である。

ma,arフィルタで生成したxがy[k]-y[k-d]差分となるようなベクトルyを生成している。

x.2 <- rnorm(10)
d <- 2
x.0 <- diffinv(x.2,d)
diff(x.0,d)

##  [1]  0.193817870  0.594084840  1.095224752  0.008422467 -0.135918666
##  [6]  2.373968878 -0.851658865  2.622396340 -1.302360586  1.481421479

x.2

##  [1]  0.193817870  0.594084840  1.095224752  0.008422467 -0.135918666
##  [6]  2.373968878 -0.851658865  2.622396340 -1.302360586  1.481421479

model$ar 引数と発散

それ以外に気になる点と言えば、arの条件として、以下のものがある。

minroots <- min(Mod(polyroot(c(1, -model$ar))))
        if (minroots <= 1) 
            stop("'ar' part of model is not stationary")

arフィルタが、累積和的な動きをすることと関係していて、再帰度が上がれば上がるほど発散に向ける力が強くなることことと対応している。

これからわかるのは、\(1 = ar[1]x+ar[2]x^2+...\)の根はすべて、その絶対値は１より大であることが、stationary(発散しない)条件であるということである。

やってみる。

この多項式条件をクリアする。 根をランダムに作成する。複素根の場合は共役複素数も根であるはず(実数係数多項式だから)であるという条件とともに、以下のように、実根数、複素根数を指定し、絶対値の１つはすべて１より大とすることにする.

これを用いてARIMAシミュレーションをする。filter()を使うバージョンと、arima.sim()を使う場合とを両方試す。

library(polynom)
n.real.roots <- 10
nhalf.complex.roots <- 5 # conjugate pairs as roots


thetas <- runif(nhalf.complex.roots)*2*pi
# All mods should be more than 1
mods <- runif(n.real.roots+nhalf.complex.roots) * 5 + 1

#ss <- sample(c(-1,1),n.real.roots+nhalf.complex.roots,replace=TRUE)

mods.real <- mods[1:n.real.roots]
mods.complex <- mods[(1+n.real.roots):(n.real.roots+nhalf.complex.roots)]
roots.R <- mods.real
roots.C1 <- mods.complex * exp(1i*thetas)
roots.C2 <- mods.complex * exp(1i*(-thetas))
roots <- c(roots.R,roots.C1,roots.C2)
polynom <- poly.from.roots(roots)

## Warning in polynomial(p): imaginary parts discarded in coercion

polynom

## 38096250000 - 159081700000*x + 311957800000*x^2 - 388603100000*x^3 +  
## 3.51642e+11*x^4 - 250031800000*x^5 + 146716300000*x^6 - 73240890000*x^7 +  
## 31698420000*x^8 - 12025670000*x^9 + 4014010000*x^10 - 1174847000*x^11 +  
## 298707500*x^12 - 65039040*x^13 + 11900750*x^14 - 1786452*x^15 +  
## 213233.1*x^16 - 19383.05*x^17 + 1255.847*x^18 - 51.48077*x^19 + x^20

ar <- (-polynom/polynom[1])[-1]
min(Mod(polyroot(c(1, -ar))))

## [1] 1.642882

x <- rnorm(1000)

x.. <- arima.sim(model=list(ar=ar),n=1000)

x. <- filter(x,ar,method="recursive")
par(mfcol=c(1,2))
plot(x.)
plot(x..)

par(mfcol=c(1,1))

rand.gen 引数

ランダム数列を発生させる際に、デフォルトでは、rnorm() 関数が使われる。

それ以外のランダム関数も使える。

x <- rexp(1000,rate=1.3)
x. <- filter(x,ar,method="recursive")

plot(x.)

x.. <- arima.sim(model=list(ar=ar),n=1000,rand.gen=rexp,rate=1.3)
plot(x..)

arima.sim()の資料でarima.sim()について確認する

http://wolfweb.unr.edu/~zal/STAT758/Lab_3.pdf　のAppendixにある例をなぞって、arima.sim()について確認する

二種類の正規乱数列Z1,Z2について、ma/convolutionタイプのフィルタリングをしている。arima.sim()ではconvolution フィルタを過去側にのみ掛けるが、リンク先資料では、両側フィルタとなっている。念のため片側の場合も付け加えてみよう。

値なしの要素が発生するが、その数に違いがあるので、シフトしてみえる。

# White noises
#===========================
Z1<-rnorm(100)
Z2<-rnorm(100,mean=10,sd=2)
# MA simulation by filter()
#===========================
X1<-filter(Z1,rep(1,10))
X2<-filter(Z2,c(.6,.3,.1))
X3 <- filter(Z1,rep(1,10),sides=1)
X4 <- filter(Z2,c(.6,.3,.1),sides=1)
par(mfcol=c(1,2))
plot(X1)
points(X3,type="l",col=2)
plot(X2)
points(X4,type="l",col=2)

par(mfcol=c(1,1))

ar/recursive タイプのフィルタリングをしてみる

# AR simulation by filter()
#===================================
X1.2<-filter(Z1,c(1,1),method='r')
X2.2<-filter(Z2,c(.6,.3,.1),method='r')
par(mfcol=c(1,2))
plot(X1.2)
plot(X2.2)

par(mfcol=c(1,1))

# ARMA simulation using arima.sim()
#===================================
par(mfrow=c(2,4))
X<-arima.sim(n=100,list(ar=c(0.9),ma=c(0.2))) # ARMA(1,1)
plot(X)
X<-arima.sim(n=100,list(ar=c(0.9))) # AR(1)
plot(X)
X<-arima.sim(n=100,list(ma=c(0.2))) # MA(1)
plot(X)
X<-arima.sim(n=100,list(ma=c(0.2)),sd=10) # MA(1) with WN(0,10)
# sd=10はrnorm()関数の引数として与えられている
plot(X)
X<-arima.sim(n=100,list(ar=c(0.9)),innov=Z1) # AR(1) with given WN Z1
# innovは乱数列(の途中から後)に特定の数列を与える。前部分はstart.innov引数で与える
plot(X)
X<-arima.sim(n=100,list(ar=c(0.9)),rand.gen=rt,df=10) # AR(1) with mild long-tail (Student's t with df=10)
plot(X)
X<-arima.sim(n=100,list(ar=c(0.9)),rand.gen=runif) # AR(1) with U[0,1] WN
plot(X)

ARMAではautocorrelation,partial autocorrelationの理論値がある

ARIMAモデルでは、タイムラグの程度ごとに、autocorrelation,partial autocorrelation の理論値がある。

# ACF and PACF computations (theoretical)
#========================================
A<-ARMAacf(ar=c(.2),ma=c(.3,.2),lag.max=20)
P<-ARMAacf(ar=c(.2),ma=c(.3,.2),lag.max=20,pacf=TRUE)
par(mfcol=c(1,2))
plot(A,type="h")
plot(P,type="h")

par(mfcol=c(1,1))

ARMAは、MAのみのモデルに変換することができる

以下では、ARMAtoMAによって、arパラメタが長さ1のベクトルであるときに、そのARMA過程がMAだったらどういうパラメタになるかを算出し、そのパラメタをma=としてARMAacfに与え、autocorrelation, partial autocorrelationが、ARMAのそれと同等になることを確認している

# MA representation for ARMA
#============================
ar <- c(0.8)
lag.max <- 20
ma.eq<-ARMAtoMA(ar=ar, lag.max=lag.max)
#P<-ARMAtoMA(ar=c(0.9,.1),ma=c(1,2,3), lag.max=5)

A<-ARMAacf(ar=ar,lag.max=lag.max)
A. <- ARMAacf(ma=ma.eq,lag.max=lag.max)
par(mfcol=c(1,2))
plot(A,A.)
abline(0,1,col=2)
P<-ARMAacf(ar=ar,lag.max=lag.max,pacf=TRUE)
P. <- ARMAacf(ma=ma.eq,lag.max=lag.max,pacf=TRUE)
plot(P,P.)
abline(0,1,col=2)

par(mfcol=c(1,1))

同じことをarパラメタベクトルの長さが2以上のときでもやってみる。

# MA representation for ARMA
#============================
ar <- c(0.3,0.1,-0.1)
lag.max <- 20
ma.eq<-ARMAtoMA(ar=ar, lag.max=lag.max)
#P<-ARMAtoMA(ar=c(0.9,.1),ma=c(1,2,3), lag.max=5)

A<-ARMAacf(ar=ar,lag.max=lag.max)
A. <- ARMAacf(ma=ma.eq,lag.max=lag.max)
par(mfcol=c(1,2))
plot(A,A.)
abline(0,1,col=2)
P<-ARMAacf(ar=ar,lag.max=lag.max,pacf=TRUE)
P. <- ARMAacf(ma=ma.eq,lag.max=lag.max,pacf=TRUE)
plot(P,P.)
abline(0,1,col=2)

par(mfcol=c(1,1))

ARIMA過程では、分散という指標があり、その理論値がある。それをmaとみなし、そのma用パラメタから算出することで近似がうまくいくことを示している。

# Variance computation
#================================
P<-ARMAtoMA(ar=0.9, lag.max=100)
sum(P^2)+1

## [1] 5.263158

ARIMAモデル推定

ARIMAもでるのar,maフィルタ推定

ARIMAがどうなっているかは分かったので、それをデータから推定する話。

arima.sim()がやっていたように、 ar部分にいくつのフィルタパラメタを置き、ma部分に、いくつを置き、時間幅として、何区分分のずれを想定するか、とによって、ARIMAデータができているものとして、推定する。

orderという引数がそれで、arパラメタ数、区分ずれ数、maパラメタ数の３整数からなる。

次の例では、order=c(1,0,0)なので、区分ずれはなく、maは関与せず、arが1変数で決まるモデルでの推定となっている。 推定は、いろいろたくさんあるが…

# ARMA estimation
#================================
X<-arima.sim(n=100,list(ar=0.9))
est<-arima(X,order=c(1,0,0))
str(est)

## List of 14
##  $ coef     : Named num [1:2] 0.87 0.86
##   ..- attr(*, "names")= chr [1:2] "ar1" "intercept"
##  $ sigma2   : num 1.15
##  $ var.coef : num [1:2, 1:2] 0.00235 0.00186 0.00186 0.60297
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2] "ar1" "intercept"
##   .. ..$ : chr [1:2] "ar1" "intercept"
##  $ mask     : logi [1:2] TRUE TRUE
##  $ loglik   : num -149
##  $ aic      : num 305
##  $ arma     : int [1:7] 1 0 0 0 1 0 0
##  $ residuals: Time-Series [1:100] from 1 to 100: 1.491 2.319 1.836 0.196 -0.17 ...
##  $ call     : language arima(x = X, order = c(1, 0, 0))
##  $ series   : chr "X"
##  $ code     : int 0
##  $ n.cond   : int 0
##  $ nobs     : int 100
##  $ model    :List of 10
##   ..$ phi  : num 0.87
##   ..$ theta: num(0) 
##   ..$ Delta: num(0) 
##   ..$ Z    : num 1
##   ..$ a    : num -1.53
##   ..$ P    : num [1, 1] 0
##   ..$ T    : num [1, 1] 0.87
##   ..$ V    : num [1, 1] 1
##   ..$ h    : num 0
##   ..$ Pn   : num [1, 1] 1
##  - attr(*, "class")= chr "Arima"

推定係数は次のように、切片(不動部分)とar1とからなっている。シミュレーションで用いたar=0.9とよく合致している

est$coef

##       ar1 intercept 
## 0.8702878 0.8601483

モデルをar+maと複雑にすると、nを大きくしないと推定値はよくならない。

# ARMA estimation
#================================
X<-arima.sim(n=100000,list(ar=c(0.2),ma=c(0.2)))
est<-arima(X,order=c(1,0,1))
est$coef

##         ar1         ma1   intercept 
##  0.21745955  0.18399060 -0.00147274

orderパラメタもいじってみよう。

# ARMA estimation
#================================
X<-arima.sim(n=100000,list(ar=c(0.2),ma=c(0.2),order=c(1,2,1)))
est<-arima(X,order=c(1,2,1))
est$coef

##       ar1       ma1 
## 0.2070198 0.1936805

autocorrelation, partial autocorrelation推定

ARIMAモデルからデータをシミュレーション作成し、その時系列データのautocorrelation, partial autocorrelationを計算し、それが、理論値とどれくらい近いかを確認してみる。

# ACF, PACF estimation
#=====================================================================
n <- 100000
ar <- c(0.3,0.2)
ma <- c(0.1)
order <- c(length(ar),0,length(ma))
lag.max <- 20
X<-arima.sim(n=n,list(ar=ar,ma=ma,order=order)) # AR(1) model
out.cor <- acf(X,lag.max = lag.max,type = c("correlation"),plot=FALSE)
theor.cor <- ARMAacf(ar=ar,lag.max=lag.max)
plot(out.cor$lag,out.cor$acf,type="h")
points(out.cor$lag,theor.cor,col=2)

out.cor2 <- pacf(X,lag.max = lag.max,type = c("correlation"),plot=FALSE)
theor.cor2 <- ARMAacf(ar=ar,lag.max=lag.max,pacf=TRUE)
plot(out.cor2$lag,out.cor2$acf,type="h")
points(out.cor2$lag,theor.cor2,col=2)