いっかくのデータサイエンティストをいく

1からプログラミングとデータサイエンスを独習したい

【時系列】GARCHモデル(R)

今回はRで実装します。 今回はこちらのサイトを写経します。

ドル円のボラティリティをGARCHで推定 – Momentum

データ

まず、データはFRED先生からドル円レートを落としてきて使います。

Japan / U.S. Foreign Exchange Rate | FRED | St. Louis Fed

Rで実装

library(tseries)
fx.ret <- fx/lag(fx)-1
fx.garch <- garch(fx.ret["1973-04::"]) #デフォルトはGARCH(1,1)
summary(fx.garch)

f:id:imakoto0323:20181204171431p:plain 結果

 
Call:
garch(x = fx.ret)

Model:
GARCH(1,1)

Residuals:
     Min       1Q   Median       3Q      Max 
-5.20882 -0.61636 -0.04313  0.55319  5.88943 

Coefficient(s):
    Estimate  Std. Error  t value Pr(>|t|)    
a0 2.120e-07   8.670e-08    2.445   0.0145 *  
a1 5.313e-02   8.385e-03    6.337 2.34e-10 ***
b1 9.421e-01   8.660e-03  108.779  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Diagnostic Tests:
    Jarque Bera Test

data:  Residuals
X-squared = 335.28, df = 2, p-value < 2.2e-16


    Box-Ljung test

data:  Squared.Residuals
X-squared = 4.8338, df = 1, p-value = 0.02791
'''

【時系列】GARCHモデル(理論編)

沖本本7章から。

はじめに

ファイナンスの世界では標準偏差のことをボラティリティとよび、重要視されています。それは分散が「最大でどれくらい損益があるか」を示していsるからです。今回はボラティリティ変動モデルです。リスクの大きさをモデリングしていきます。

ARCHモデル

ARCHモデルは自己回帰条件付分散不均一モデルのことです。 このモデルの考え方は「絶対値の大きなノイズが前回来たら、今回の分散は大きくなる」というものです。

1次のARCHモデルは以下のように定式化できます。

f:id:imakoto0323:20181204163205p:plain

日本語訳すると データ=期待値+ノイズ ノイズ=√条件付き分散×分散1のホワイトノイズ 条件付き分散=ω+a1×(前期のノイズ)の二乗

となります*1

GARCHモデル

ARCHモデルをより一般化したモデルがGARCHモデルです。訳すと一般化自己回帰条件付き分散不均一モデルとなります。 このモデルでは「より長くデータのブレが広がる状況が持続するモデルを少ないパラメータで表現したい」という目的があります。

GARCH(1,1)の定式は

f:id:imakoto0323:20181204164019p:plain

となります。 これも書き下すと データ=期待値+ノイズ ノイズ=√条件付き分散×分散1のホワイトノイズ 条件付き分散=ω+α1×(前期のノイズ)の二乗+β1×前期の条件付き分散

変わったのは3番目の式だけです*2

*1:なぜこの式で表現できているのかは分かりません

*2:これも正直よく分かりません。誰か教えてください笑

【時系列】共和分

沖本本5章から。

共和分

単位根の持つデータ同士で回帰分析した場合、見せかけの回帰になってしまうことが多いです。もうあきらめるしかないのか。。。 そんなことはありません!見せかけではない場合があります! それは共和分を持っているときです。 共和分の概念自体は単純でf:id:imakoto0323:20181204131239p:plainf:id:imakoto0323:20181204131311p:plainが線形結合すなわち

f:id:imakoto0323:20181204131404p:plain

の場合に単位根を持たなくなったとき、両者は共和分の関係にあるといいます。共和分の関係の場合、回帰分析ができます

共和分検定

共和分の検定にはEngle-Granger検定があります。 この検定の考え方はシンプルで ①OLSで回帰分析 ②残差が単位根がなければ共和分ありとみなす この検定の帰無仮説は「共和分関係がない」です。

Rで実装

#パッケージ
library(vars)


#ランダムウォークを生成
n_sample <- 400 #サンプルサイズ
set.seed(10)
rw <- cumsum(rnorm(n=n_sample))
x_co <- 0.6*rw + rnorm(n=n_sample)
y_co <- 0.4*rw + rnorm(n=n_sample)

#ADF(単位根)検定
summary(ur.df(y_co, type = "none" ))
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
  lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
  Min      1Q  Median      3Q     Max 
-4.0281 -0.8282 -0.0794  0.8101  3.4312 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -0.012165   0.009448  -1.288    0.199    
z.diff.lag -0.453423   0.044548 -10.178   <2e-16 ***
  ---
  Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.236 on 396 degrees of freedom
Multiple R-squared:  0.2162, Adjusted R-squared:  0.2122 
F-statistic: 54.62 on 2 and 396 DF,  p-value: < 2.2e-16


Value of test-statistic is: -1.2875 
#単位根を持つという帰無仮説は棄却できない
Critical values for test statistics: 
  1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

summary(ur.df(x_co, type = "none" ))
############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
  lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
  Min      1Q  Median      3Q     Max 
-4.6304 -1.0182 -0.0304  0.8609  4.3365 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
z.lag.1    -0.007435   0.007294  -1.019    0.309    
z.diff.lag -0.418492   0.045650  -9.167   <2e-16 ***
  ---
  Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.436 on 396 degrees of freedom
Multiple R-squared:  0.1803, Adjusted R-squared:  0.1762 
F-statistic: 43.57 on 2 and 396 DF,  p-value: < 2.2e-16


Value of test-statistic is: -1.0194 
#単位根を持つという帰無仮説は棄却できない
Critical values for test statistics: 
  1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62

#データを1つにまとめる
df <- data.frame(
  y_co = y_co,
  x_co = x_co,
)

#ts型に変換
ts_df <- ts(df)

#共和分検定
summary(ca.po(ts_df, demean = "none"))
######################################## 
# Phillips and Ouliaris Unit Root Test # 
######################################## 

Test of type Pu 
detrending of series none 


Call:
  lm(formula = z[, 1] ~ z[, -1] - 1)

Residuals:
  Min         1Q     Median         3Q        Max 
-2.580e-15 -3.550e-16  4.000e-18  3.650e-16  5.676e-14 

Coefficients:
  Estimate Std. Error    t value Pr(>|t|)    
z[, -1]x_co  6.667e-01  1.593e-17  4.184e+16   <2e-16 ***
  z[, -1]z    -6.667e-01  8.722e-17 -7.644e+15   <2e-16 ***
  ---
  Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.133e-15 on 398 degrees of freedom
Multiple R-squared:      1,  Adjusted R-squared:      1 
F-statistic: 8.822e+32 on 2 and 398 DF,  p-value: < 2.2e-16


Value of test-statistic is: -1.357322e+16 
#帰無仮説を棄却できない
Critical values of Pu are:
  10pct    5pct    1pct
critical values 26.7022 32.9392 46.4097

これまでのVAR~単位根までの過程についてこちらに詳しくまとめられています。

tjo.hatenablog.com

【時系列】差分系列と単位根検定

沖本本5章から。

差分系列

非定常なデータは時間によって期待値が変化するため予測ができません。 しかし、非定常なデータでも差分をとることによって定常になることがあります。これを差分系列と呼びます。差分をとることはよくやることで一見分析ができないものでも分析ができるようになります。 f:id:imakoto0323:20181204124025p:plain

単位根と見せかけの回帰

差分をとった場合、原系列では非定常だったデータが定常になる場合があります。これを単位根過程といいます。

単位根過程で何が困るのかというと、時系列データの回帰分析をする場合です。時系列データを回帰する場合、本来は因果関係のないデータ同士でも因果関係があるかのように検出されます。これを見せかけの回帰といいます。 なぜそうなるかというと、上昇、下降が1/2の確率で発生すると(上昇、上昇)(下降、下降)が何回か重なって発生するためかなと個人的には思っています。

Rで実装

library(tseries)
x <- rnorm(1000)  # サンプルとして定常過程を使う
adf.test(x)
> Augmented Dickey-Fuller Test

data:  x
Dickey-Fuller = -8.7153, Lag order = 9, p-value = 0.01
alternative hypothesis: stationary
# 帰無仮説が棄却された、すなわちこれは定常過程という判定


y <- diffinv(x) # 単位根を含むように変えてみた
adf.test(y)

Augmented Dickey-Fuller Test

data:  y
Dickey-Fuller = -3.121, Lag order = 9, p-value = 0.1038
alternative hypothesis: stationary
# 帰無仮説が採択された、すなわちこれは単位根過程

【時系列】ホワイトノイズとランダムウォークの違いを考えてみる

ふとホワイトノイズのランダムウォークってどう違うんだろう?と考えてみたのでこの2つをまとめてみます

ホワイトノイズ

ホワイトノイズは

f:id:imakoto0323:20181204113547p:plain となり弱定常性となります。

ランダムウォーク

時刻nにおける価格p(n)を、

p(n)=p(n−1)+d(n)

のように1サンプル前の価格p(n−1)にd(n)(ただし、d(n)は1か−1)を加えて算出していく時系列を単純ランダムウォークといいます。ランダムウォークは定常性を持ちません。

両者の違いは

一定のところでランダムに動くのがホワイトノイズでどこへ向かうのかわからないのがランダムウォークというところでしょうか。

【時系列】VAR(インパルス応答)

沖本本4章3から。

インパルス応答とは

ある変数にインパクトを与えると、その影響がどれくらい続くのかをインパルス応答といいます。消費が急に増えると収入にどんな影響があるのかを定量的に評価できます。 そして「時間遅れ(タイムラグ)」の向きを見ることで、「どの時系列がどの別の時系列に対して影響を与える=どちらが原因でどちらが結果か」が分かります。

Rで実装

irf関数でサクッと求められます。

msci.irf<-irf(Canada.var,n.ahead=14,ci=0.95)
# 14期先までインパルス応答関数を計算
# 信頼区間はブートストラップ法で算出される
plot(msci.irf)

f:id:imakoto0323:20181204105430p:plain

【時系列】VAR(グレンジャーの因果性検定)

沖本本4章3から。

グランジャーの因果性

現在と過去のxの値だけに基づいた将来のxの予測と、現在と過去のxとyの値に基づいた将来のxの予測を比較して、後者のMSE(残渣平方和)の方が小さくなる場合、ytからxtへのグレンジャー因果性(Granger causality)が存在するといわれています。

つまり、個人収入と個人消費の場合、

今年の個人収入 = 去年の個人収入×傾きA + 去年の個人消費×傾きB + 切片C

今年の個人収入 = 去年の個人収入×傾きA + 切片C

を比較した場合、前者のほうがMSEが小さい場合、去年の個人消費を変数に入れたほうがいいということになります。 f:id:imakoto0323:20181031211009p:plain

Granger因果による 時系列データの因果推定(因果フェス2015)から拝借)

グランジャーの因果性検定

グランジャーの因果性検定とは、現在と過去のxの値だけに基づいた将来のxの予測と、現在と過去のxとyの値に基づいた将来のxの予測を比較して、後者のMSE(残渣平方和)の方が有意に小さくなるか検定するものです。

グランジャーの因果性検定の手順は (1) VARモデルにおけるyktのモデルをOLSで推定し、その残差平方和をSSR1とする。

(2) VARモデルにおけるyktのモデルに制約を課したモデル(例でいうと個人消費がないモデル)をOLSで推定し、その残差平方和をSSR0とする。

(3) F統計量を F≡(SSR0−SSR1)/rSSR1/(T−np−1) で計算する。ここで、rはグレンジャー因果性検定に必要な制約の数である。

(4) rFをχ2(r)の95%点と比較し、rFのほうが大きければ、ある変数(群)からyktへのグレンジャー因果性は存在し、小さければグレンジャー因果性は存在しないと結論する。

という方法です。

Rで実装

Rで計算するにはcausality(){vars}関数を使います。

causality(Canada.var, cause = "e")

$`Granger`

    Granger causality H0: e do not Granger-cause prod rw
    U

data:  VAR object Canada.var
F-Test = 5.5053, df1 = 9, df2 = 240, p-value =
6.882e-07


$Instant

    H0: No instantaneous causality between: e and prod
    rw U

data:  VAR object Canada.var
Chi-squared = 23.585, df = 3, p-value = 3.049e-05