0.0 0.1 0.2 0.3 0.4 0.5
0.00.20.40.60.81.0
frequency Series: lh
0.0 0.1 0.2 0.3 0.4 0.5
0.00.20.40.60.81.0
frequency AR(3) fit to lh
0 1 2 3 4 5 6
0.00.20.40.60.81.0
frequency Series: ldeaths
黄体形成ホルモン量データの累積 ペリオドグラム
原データ(上左)
AR(3) モデルを当てはめた残差
(上右)
英国の肺疾患死亡者数(両性)デー タの累積ペリオドグラ(下左)
demean = TRUE, series, ...)
#
クラス"ar"
に対するS3
メソッドpredict(object, newdata, n.ahead = 1, se.fit = TRUE, ...)
引数：
x
一変量，もしくは多変量時系列aic
論理値フラグ．もしTRUE
なら，自己回帰モデルの次数を選択するために赤池の 情報量規準が使われる．もしFALSE
なら次数order.max
のモデルが当てはめら れるorder.max
当てはめられるモデルの(
最大)
次数．既定ではN
を観測数としてN-1
と10*log10(N)
の小さいほうが使われる．method="mle"
では，この数と12
の小 さいほうとされるmethod
モデルの当てはめに使われる手法を与える文字列．既定の引数のどれかである必要があるが，最初の数文字だけで十分である．既定では
"yule-walker"
na.action
欠損値を処理するために使われる関数名demean
論理値．当てはめの過程で平均値を推定するか?
series
時系列の名前．既定ではdeparse(substitute(x))
var.method
イノベーション分散を推定する手法...
特定の手法のためのオプション追加引数object ar()
関数による当てはめ結果newdata
予測を試みるデータn.ahead
何ステップ先まで予測を行うか?
se.fit
論理値．予測誤差の標準誤差の推定値を返すか?
返り値： 以下の成分を持つクラス
"ar"
のリスト：order
当てはめられたモデルの次数．aic=TRUE
ならば，AIC
を最小化するものが，さもなければ
order.max
そのものar
当てはめられたモデルの自己回帰係数の推定値var.pred
予測分散．時系列の分散のうち，自己回帰モデルで説明できなかった部分の推定値
x.mean
時系列の平均値の推定値で，当てはめ過程と予測で用いられるx.intercept ar.ols()
専用．x-x.mean
に対するモデルの切片値aic aic
引数の値n.used
時系列中の総観測値数order.max order.max
引数の値partialacf
最大ラグorder.max
までの部分自己相関系数の推定値resid
「次数1
」の観測値で条件付けられた，当てはめモデルの残差．「次数1
」の残差は
NA
とされる．もしx
が時系列ならresid
も時系列になるmethod
引数method
の値series
時系列の名前frequency
時系列の周波数call
マッチした呼出し式asy.var.coef (
一次元時系列で，order>0
の場合)
漸近理論に基づく系数推定値の分散行列
predict.ar()
，予測値時系列，もしくはse.fit=TRUE
に対しては，予測値である成分pred
と推定標準誤差である成分se
を持つリストである．成分はともに時系列になる．時系列
{ x
t}
に対するp
次のAR
モデルAR(p)
は次の式で定義される: (x
t− m) = a
1(x
t−1− m) + . . . + a
p(x
t−p− m) + ε
tここで，
m
は平均値，a
1, . . . , a
p は自己回帰係数，そしてε
tは時点t
での誤差項である．関 数
ar()
は 様 々 な 手 法 で 自 己 回 帰 モ デ ル の 当 て は め を 行 う 関 数 群 ，ar.yw(), ar.burg(), ar.ols()
そしてar.mle()
，へのラッパ関数である．これらは，aic=TRUE
な らばAIC
法でモデルの選択を行う．しかしながら，真正の最尤推定法を用いるar.mle()
を除けばこれは問題なしとしない．AIC
値は分散の推定値が最尤推定値であるかのよう に計算され，尤度から行列式項は除かれる．これは，推定パラメータを用いて評価された 正規分布尤度とは異なるものであることを注意しよう．ar.yw()
では，イノベーション は当てはめ係数とx
の自己共分散から計算される．ar.burg()
ではイノベーション分散，したがってAIC
値，の計算に二つの手法を用いることができる．手法
1
はS-PLUS
と同様にLevinson-Durbin
の更新式を更新のために 用いる．手法2
は，前進・後退予測誤差の2
乗和の平均を用いる．推定係数は分散の推定 法に(
多少)
依存するであろう．ar()
関数は，x
の全体平均を引き去るか，引き去る定数を推定(ar.mle())
すること により，既定でモデル中に定数を含むことを注意しよう．注意：
ar.mle()
は一次元時系列に対してだけ実装されている．method="mle"
のよる 長い時系列当てはめは非常に時間がかかる可能性がある．関連：
ar.ols(), ARMA
モデルに対してはarima0()
，自己相関係数からのAR
モデル の構築にはacf2AR()
．# 黄体形成ホルモン量データlhを使用
> ar(lh) # 既定のYule-Walker法ar.yw()による当てはめ
Call:
ar(x = lh) # 呼出し式
Coefficients:
1 2 3 # 当てはめモデルの係数
0.6534 -0.0636 -0.2269
Order selected 3 sigma^2 estimated as 0.1959 # 3次のARモデル選択
> str(ar(lh)) # 当てはめ結果の詳細
List of 14
$ order : int 3
$ ar : num [1:3] 0.6534 -0.0636 -0.2269
$ var.pred : num 0.196
$ x.mean : num 2.4
$ aic : Named num [1:17] 18.307 0.996 0.538 0.000 1.490 ...
..- attr(*, "names")= chr [1:17] "0" "1" "2" "3" ...
$ n.used : int 48
$ order.max : num 16
$ partialacf : num [1:16, 1, 1] 0.576 -0.223 -0.227 0.103 -0.076 ...
$ resid : Time-Series [1:48] from 1 to 48: NA NA NA -0.200 -0.169 ...
$ method : chr "Yule-Walker"
$ series : chr "lh"
$ frequency : num 1
$ call : language ar(x = lh)
$ asy.var.coef: num [1:3, 1:3] 0.02156 -0.01518 0.00482 -0.01518 0.03117 ...
- attr(*, "class")= chr "ar"
> ar(lh, method="burg") # ar.burg()使用 Call:
ar(x = lh, method = "burg") Coefficients:
1 2 3
0.6588 -0.0608 -0.2234
Order selected 3 sigma^2 estimated as 0.1786
> ar(lh, method="ols") # ar.ols()使用 Call:
ar(x = lh, method = "ols") Coefficients:
1 0.586
Intercept: 0.006234 (0.06551)
Order selected 1 sigma^2 estimated as 0.2016
> ar(lh, FALSE, 4) # 4次のARモデル当てはめ
Call:
ar(x = lh, aic = FALSE, order.max = 4) Coefficients:
1 2 3 4
0.6767 -0.0571 -0.2941 0.1028
Order selected 4 sigma^2 estimated as 0.1983
> (sunspot.ar <- ar(sunspot.year)) # 太陽黒点数データsunspot使用 Call:
ar(x = sunspot.year) Coefficients:
1 2 3 4 5 6 7 8
1.1305 -0.3524 -0.1745 0.1403 -0.1358 0.0963 -0.0556 0.0076 9
0.1941
Order selected 9 sigma^2 estimated as 267.5 # 9次のARモデル選択
> predict(sunspot.ar, n.ahead=25) # 当てはめモデルによる予測
$pred # 予測値時系列
Time Series:
Start = 1989 End = 2013 Frequency = 1
[1] 135.25933 148.09051 133.98476 106.61344 71.21921 40.84057 18.70100 [8] 11.52416 27.24208 56.99888 87.86705 107.62926 111.05437 98.05484 [15] 74.84085 48.80128 27.65441 18.15075 23.15355 40.04723 61.95906 [22] 80.79092 90.11420 87.44131 74.42284
$se # 予測標準誤差時系列
Time Series:
Start = 1989 End = 2013 Frequency = 1
[1] 16.35519 24.68467 28.95653 29.97401 30.07714 30.15629 30.35971 30.58793 [9] 30.71100 30.74276 31.42565 32.96467 34.48910 35.33601 35.51890 35.52034 [17] 35.65505 35.90628 36.07084 36.08139 36.16818 36.56324 37.16527 37.64820 [25] 37.83954
# 2 変量時系列への当てはめ例．売上高データBJsales使用
> ar(cbind(BJsales, BJsales.lead), method="burg") Call:
ar(x = cbind(BJsales, BJsales.lead), method = "burg")
$ar # 4次の2変量ARモデルを選択
, , 1 # 以下各次数の(行列)係数
BJsales BJsales.lead BJsales 1.21197 0.07312 BJsales.lead 0.07411 0.45684 , , 2
BJsales BJsales.lead
BJsales -0.26022 -0.1120
BJsales.lead -0.06904 0.3111 , , 3
BJsales BJsales.lead
BJsales -0.01754 3.93591 BJsales.lead 0.01792 0.09632 , , 4
BJsales BJsales.lead BJsales -0.07746 -1.33836 BJsales.lead -0.01158 -0.09118
$var.pred
BJsales BJsales.lead BJsales 0.38426 0.01315 BJsales.lead 0.01315 0.07657
> ar(BJsales, method="burg") # 1変量で当てはめるとAR(5)を選択 Call:
ar(x = BJsales, method = "burg") Coefficients:
1 2 3 4 5
1.2213 -0.0619 -0.0649 0.0668 -0.1668 Order selected 5 sigma^2 estimated as 1.759
4.4.2 最小自乗法による AR モデルの当てはめ ar.ols()
通常最小自乗法
(OLS, Ordinary Least Square)
により，時系列にAR
モデルを当て はめる．既定ではAIC
法で次数を決める．書式：
ar.ols(x, aic = TRUE, order.max = NULL, na.action = na.fail, demean = TRUE, intercept = demean, series, ...)
引数：
x
一(
多)
変量時系列aic
論理値フラグ．もしTRUE
なら，AIC
法を用いてAR
モデルの次数を選択する．もし
FALSE
なら次数order.max
のモデルを使用するorder.max
当 て は め モ デ ル の(
最 大)
次 数 ．既 定 で は ，観 測 値 数 をN
と す る と10*log10(N)
が使われるna.action
欠損値を処理する関数demean x
から平均値を引き去ったAR
モデルを当てはめるべきか？intercept
切片(
定数)
項を別個に当てはめるべきか？series
時系列の名前．既定ではdeparse(substitute(x)) ...
他のメソッドへ(
から)
引き渡される追加引き数返り値： クラス
"ar"
オブジェクトで，以下の成分を持つ：order
当てはめモデルの次数．これは，もしaic=TRUE
ならAIC
の最小化で選ばれる．さもなければ
order.max
であるar
当てはめモデルのAR
係数var.pred
予測分散．つまり，時系列の分散のうち，AR
モデルで説明できない部分x.mean
当てはめで使われた推定平均(
もしdemean=FALSE
なら0)
で，予測でも使わ れるx.intercept x-x.mean
に対するモデルの切片項．intercept=FALSE
なら0 aic aic
引き数の値n.used
時系列中の観測値数order.max order.max
引き数値partialacf NULL
．ar()
との互換性のためにあるresid
第「一次」の観測値で条件付けた，当てはめモデルからの残差．第「一次」の残差は
NA
とされる．もしx
が時系列なら，rsid
も時系列method
文字列"Unconstrained LS"
series
時系列の名前asy.se.coef
係数推定値の漸近理論による標準誤差ar.ols()
は一般のAR
モデルを時系列x
に当てはめる．多変量時系列でも，可能性として非定常でも良い．たとえ，時系列の一部分が非定常でも，また
(
同じ階数で)
階差 をとってあっても，結果の制約無し最小自乗推定値は一致性を持つ．正確にはAR
係数 は次のように符号が決められている：(x
t− m) = a
0+ a
1(x
t−1− m) + . . . + a
p(x
t−p− m) + e
tここで
a
0 はintercept = TRUE
でない限り0
であり，m
はdemean = TRUE
なら標本 平均，さもなければ0
である．次数の選択は，もしaic = TRUE
ならAIC
法で行われ る．これは，ar.ols()
が真の最尤推定を行わないので，問題無しとはしない．AIC
は(
残差の分散行列から計算された)
分散推定値があたかもMLE
であるかのように計算さ れ，尤度から行列式項を除外してある．これは，推定パラメータ値を用いて計算した正規 分布尤度とは異なるものであることを注意しよう．もしintercept = TRUE
でdemean
= FALSE
であるときは，少し注意がいる．時系列がほぼ0
の周りに中心化されているときだけこれを使用しよう．さもなければ，計算は不正確になるか，全く失敗するかも知れ ない．
# 黄体形成ホルモン量データlhを使用(以下の図を参照)
> (ts <- ar(lh, method="burg")) # ar()法によるARモデルの当てはめ
> ts.plot(lh, predict(ts, n.ahead = 10)$pred) # 時系列と予測値の同時プロット
> (ts <- ar.ols(lh)) # ar.ols()法でARモデル当てはめ
Call:
ar.ols(x = lh) Coefficients:
1 0.586
Intercept: 0.006234 (0.06551)
Order selected 1 sigma^2 estimated as 0.2016
> ts.plot(lh, predict(ts, n.ahead = 10)$pred) # 時系列と予測値の同時プロット
> ts <- ar.ols(lh, FALSE, 4) # AR(4)の当てはめ
> ts.plot(lh, predict(ts, n.ahead = 10)$pred) # 時系列と予測値の同時プロット
> ts <- ar.ols(lh, order.max=6, demean=FALSE, intercept=TRUE)
> data(BJsales) # 企業売り上げデータ読み込み
> ts <- ar.ols(BJsales) # ARモデルのの当てはめ
> ts.plot(BJsales, predict(ts, n.ahead=10)$pred) # 系列，予測値の同時プロット
# ヨーロッパ株価指標データEuStockMarketsを使用
> x <- diff(log(EuStockMarkets)) # 対数値の一階階差を取る
# 最大6次までのARモデルの当てはめ，平均化せず，定数項を含める
> ar.ols(x, order.max=6, demean=FALSE, intercept=TRUE) Call:
ar.ols(x = x, order.max = 6, demean = FALSE, intercept = TRUE)
$ar , , 1
DAX SMI CAC FTSE
DAX 0.004560 -0.095781 0.039975 0.04856 SMI -0.009204 -0.007142 0.037758 0.06826 CAC -0.026624 -0.113688 0.063807 0.09154
FTSE -0.010299 -0.089246 -0.003195 0.16409
$x.intercept
DAX SMI CAC FTSE
0.0006941 0.0007813 0.0004866 0.0004388
$var.pred
DAX SMI CAC FTSE
DAX 1.056e-04 6.683e-05 8.274e-05 5.192e-05 SMI 6.683e-05 8.496e-05 6.252e-05 4.254e-05 CAC 8.274e-05 6.252e-05 1.207e-04 5.615e-05 FTSE 5.192e-05 4.254e-05 5.615e-05 6.224e-05
Time
0 10 20 30 40 50 60
1.52.02.53.03.5
Time
0 10 20 30 40 50 60
1.52.02.53.03.5
Time
0 10 20 30 40 50 60
1.52.02.53.03.5
Time
0 50 100 150
200210220230240250260
黄 体 形 成 ホ ル モ ン 量 デ ー タ へ の ARモデル当てはめと 10期先ま での予測
(上左)ar()による当てはめ (上右)ar.ols()による当てはめ (下左) AR(4)の当てはめ (下右)ar.ols() による企業売り 上げデータへのARモデルの当て はめ10期先までの予測