状態空間モデルとカルマンフィルタ .1 時系列の構造モデル StructTS().1時系列の構造モデルStructTS()

時系列に構造モデル

(structual model)

を当てはめる．

書式：

StructTS(x, type = c("level", "trend", "BSM"), init = NULL, fixed = NULL, optim.control = NULL)

引数：

x

一変量数値時系列．欠損値があっても良い

type

構造モデルのクラス．もし与えられないと，

BSM

が

frequency(x)>1

の時系 列に，さもなければ，局所傾向モデルが当てはめられる

init

分散パラメータの初期値

fixed

パラメータの総数と同じ長さのオプションの数値ベクトル．もし与えられると，

fixed

中の

NA

でない項目だけが変化する．おそらく最も有用な使い途は，分散

を

0

と設定すること

optim.control optim()

に対する制御パラメータのリスト．

"L-BFGS-B"

法使用 返り値： クラス

"StructTS"

のリストで，次の成分を持つ：

coef

成分の推定分散値

loglik

最大対数尤度値．これらの全てのモデルは非定常であり，これは幾つかの観測

値に対する拡散プライア

(diffuse prior)

を含み，従って

arima()

や異なったタ イプの構造モデルと比較はできない

data

時系列

x

residuals

標準化残差

fitted

多変量時系列で，水準，傾き，季節成分に対する成分を持つ．これらは経時的

(

時系列の最後ではなく，時間

t

毎に

)

に推定される

call

呼出し式

series

時系列

x

の名前

convergence optim()

が返した値

model,model0

当てはめで使われたカルマンフィルタを表すリスト．

KalmanLike()

を見よ．

model0

がフィルタの初期値であり，

model

が最終状態

xtsp x

の

tsp

属性

構造時系列モデルは，時系列の幾つかの成分への分解に基づく，一次元時系列の

(

線形 ガウシアン

)

状態空間モデルである．これは，幾つかの誤差分散により指定され，そのう ちの幾つかは

0

であっても良い．

最も簡単なモデルは局所水準モデルで，

type = "level"

により指定される．これは 次式で展開する背景水準

m

t を持つ：

m

t+1

= m

t

+ ξ

t

, ξ

t

∼ N(0, σ

²

ξ).

観測値のモデルは次式で与えられる：

x

t

= m

t

+ ε

t

, ε

t

∼ N(0, σ

²

ε).

二つのパラメータ

σ

²

ξ

と

σ

²

ε

がある．これは

ARIMA(0,1,1)

モデルであるが，パラメー タセットに制約がある．

局所線形傾向モデル

type = "trend"

は同じ測定方程式を持つが，

m

t 系列に次式

(3

つの分散パラメータを持つ

)

のように時間とともに変化する傾きを持つ：

m

t+1

= m

t

+ n

t

+ ξ

t

, ξ

t

∼ N(0, σ

_ξ²

), n

t+1

= n

t

+ ζ

t

, ζ

t

∼ N(0, σ

²_ζ

).

σ

_ζ²

= 0 (

局所水準モデルに帰着

)

や滑らかな傾向を意味する

σ

²_ξ

= 0

に成ることも珍しく ない．これは制約された

ARIMA(0,2,2)

モデルである．

基本の構造モデル

type = "BSM"

は追加の季節成分を持つ局所傾向モデルである．

従って計測方程式は

x

t

= m

t

+ s

t

+ ε

t

, ε

t

∼ N(0, σ

_ε²

)

となり，ここで

s

tは次の式に従い変化する季節成分である：

s

t+1

= − s

t

− . . . − s

t−s+2

+ w

t

, w

t

∼ N(0, σ

²_w

).

境界ケース

σ

²_w

= 0

は非ランダム

(

しかし任意の

)

季節パターンに相当する．

(

これはしば しば

BSM

の「ダミー変数」版と呼ばれる．

)

注意：構造モデルの当てはめは，多くの参考文献がそう述べているよりはるかに困難であ る．例えば，

Brockwell & Davis

は

AirPassengers

データを議論しているが，彼らの解 は局所的最大解のように思われ，

StructTS()

が与える当てはめ程は良くは無い．一つも しくは複数の変数が

0

である当てはめを見出すことはごく普通であり，これは

σ

_ε²も例外 ではない．

# 木の年輪幅データtreeringを使用(以下の図を参照)

> trees <- window(treering, start=0) # 紀元後のデータだけ取り出す

> fit1 <- StructTS(trees, type = "level") # 構造モデル当てはめ

> tsdiag(fit1) # 年輪データ当てはめの診断図

> plot(cbind(fitted(fit1), resids=resid(fit1)), main = "treering")

# 英国のガス消費量データUKgasを使用(以下の図を参照)

> fit2 <- StructTS(log10(UKgas), type = "BSM") # 構造モデル当てはめ

> plot(cbind(fitted(fit2), resids=resid(fit2)), main = "UK gas consumption")

0.80.91.01.11.21.3

fitted(fit1) −4−3−2−1012

resids

0 500 1000 1500 2000

Time

UK gas consumption Standardized Residuals

Time

0 500 1000 1500 2000

−4−3−2−1012

0 5 10 15 20 25 30

0.00.20.40.60.81.0

Lag

ACF

ACF of Residuals

2 4 6 8 10

0.00.20.40.60.81.0

p values for Ljung−Box statistic

lag

p value

(左)年輪データへの構造モデル当てはめ結果．上から水準，残差系列．(右)年輪データ への構造モデル当てはめ結果の診断図

2.02.22.42.62.8

fitted(fit2).level −0.04−0.020.000.02

fitted(fit2).slope −0.4−0.20.00.10.20.3

fitted(fit2).sea −2024

resids

1960 1965 1970 1975 1980 1985

Time UK gas consumption

(下左) 英国ガス消費量への構造モデルの当 てはめ結果．上から，水準，傾き，周期，残 差系列

4.11.2 固定区間平滑化 tsSmooth()

状態空間モデルを用い，一変量時系列に固定区間平滑化を実行する．固定区間平滑化 は，全ての観測値に基づき，各時点での最良の状態の推定を与える．

書式：

tsSmooth(object, ...)

引数：

object

当てはめ時系列．現在クラス

"StructTS"

のオブジェクトだけがサポートさ れている

...

将来のメソッドへの可能な引き数

返り値： 一次元時系列で，原時系列の状態空間と各時点での結果に等しい次元を持つ．

(

季節モデルに対しては，現在の季節成分だけが返される

)

# 国際航空路線総顧客数データAirPassengersを使用(以下の図を参照)

> ap <- log10(AirPassengers) - 2

> fit1 <- StructTS(ap, type= "BSM") # 構造モデルの当てはめ

> par(mfrow=c(1,2))

# 原データと構造モデルによる当てはめ系列のプロット

> plot(cbind(ap, fitted(fit1)), plot.type = "single")

# 原データと固定区間平滑化のプロット

> plot(cbind(ap, tsSmooth(fit1)), plot.type = "single")

# 企業の利益データJohnsonJohnsonを使用(以下の図を参照)

> JJ <- log10(JohnsonJohnson)

> plot(JJ) # 原データのプロット

> fit2 <- StructTS(JJ, type="BSM") # 構造モデルの当てはめ

> sm <- tsSmooth(fit2) # 固定区間平滑化

# 水準，傾き，季節変動のプロット

> plot(cbind(sm[, 1], sm[, 2], sm[, 3]), plot.type = "single", col = c("black", "green", "blue"))

> abline(h = -0.5, col = "grey60")

Time

cbind(ap, fitted(fit))

1950 1954 1958

0.00.20.40.60.8

Time

cbind(ap, tsSmooth(fit))

1950 1954 1958

0.00.20.40.60.8 JJ

19601965 197019751980

0.00.51.0 cbind(sm[, 1], sm[, 2], sm[, 3])

19601965197019751980

−0.20.00.20.40.60.81.01.2

固定区間平滑化．航空路線顧客数データ．

(左) 原データ，構造モデルと残差．(右)原 データ，固定区間平滑化と残差

固定区間平滑化．企業売上データ．(左) 原 データ．(右)固定区間平滑化による水準，傾 向，季節系列

4.11.3 構造モデル用補助関数 KalmanLike()

カルマンフィルタを用いて，

(

ガウシアン

)

対数尤度を求める，また予測や平滑化を行う．

書式：

KalmanLike(y, mod, nit = 0) KalmanRun(y, mod, nit = 0) KalmanSmooth(y, mod, nit = 0) KalmanForecast(n.ahead = 10, mod)

makeARIMA(phi, theta, Delta, kappa = 1e6)

引数：

y

一変量時系列

mod

状態空間モデルを記述するリスト．以下の説明を参照

nit

初期化が計算される時間．

nit=0

は初期化が一時点先の予測様であることを意 味し，したがって

Pn

は最初のステップでは計算すべきではないことを意味する

n.ahead

予測が必要な将来のステップ数

phi,theta 0

以上の長さの数値ベクトルで，

AR

と

MA

パラメータを与える

Delta

階差係数のベクトルで，

ARMA

モデルが

y[t]-Delta[1]*y[t-1]-...

に当 てはめられる

kappa

階差モデル中の過去の観測値に対する事前分散

(

イノベーション分散への割合

で与える

)

これらの関数は，状態ベクトル

a

，推移

a ← T a + e, e ∼ N(0, κQ)

，観測値式

y = Z

^′

a + Rη, η ∼ N(0, κh).

を持つ一般の一変量状態空間モデルに関連する．尤度は

κ

を推定した後で得られるプロファイル尤度^*6 である．モデルは少なくとも次の成分を持 つリストで指定される：

T

推移行列

Z

観測値係数

h

観測値分散

V RQR

a

現在の状態推定

P

状態の不確定性行列の現在の推定値

Pn

状態の不確定性行列の時間

t-1

の推定値

KalmanSmooth

クラス

"tsSmooth"

に対する実動関数

makeARIMA ARIMA

モデルに対する状態空間モデルを構成する

各関数の返り値リストの成分は以下のようになる：

KalmanLike

成分

Lik(

尤度からある定数を引いたもの

)

と，

κ

の推定値である

s2 KalmanRun KalmanLike()

の出力である長さ

2

のベクトル

values

，残差

resid

，各

時点に対して一行を持つ行列である経時的状態推定値である

states

KalmanSmooth

成分

smooth

は全ての観測値に基づく状態推定の

nxp

行列で，各時点 に対して一つの行を持つ．成分

var

は分散行列から成る

nxpxp

配列

KalmanForecast

予測値を与える成分

pred

と

( s2

倍された

)

スケール化されていな い予測誤差

makeARIMA

その引き数を成分として持つモデルリスト

注意：これらの関数は，引き渡された引き数の整合性をチェックする他の関数から呼び出 されるようにデザインされており，自分自身ではほとんどチェックを行わない．

ドキュメント内 ためになった他の人のサイト script of (ページ 143-147)