2017-08-31

日本国債金利の主成分分析

概要

今更ながら金利モデルに興味が出てきたので，記事を書いてみた
主成分分析を用いて，日本国債金利の構造を推計してみた
（もちろん，今回も内容の信ぴょう性はありません）

Nelson-Siegelモデルについて

　金利の期間構造に関するモデルは，たくさんある *1

　Nelson-Siegelモデル*2は，その代表的なモデルである。このモデルは，スポットレートカーブの構造を比較的に単純に説明することができる。

　残存期間mのフォワードレートをr(m)とする。

$r(m) = \beta_{0} + \beta_{1} \exp( - m / \tau ) + \beta_{2} [ (m / \tau ) \exp ( - m / \tau ) ]$

　残存期間mのスポットレートをR(m)とすると，R(m)はフォワードレートの平均値である。

$R(m) = 1/m \int_{0}^{m} r(x) dx$

　この2つの式を組み合わせる。

$R(m) = \beta_{0} + \beta_{1} [ \frac{1 - \exp ( -m / \tau )}{ ( m / \tau ) } ] + \beta_{2} [ \frac{ 1 - \exp ( - m / \tau ) }{ (m / \tau ) } - \exp( - m / \tau ) ]$

ここで，各パラメータβ0，β1，β2，τが時変的（time-varying）であると仮定すると，ダイナミックなNelson-Siegelモデルとなる。

$R(m) = \beta_{0t} + \beta_{1t} [ \frac{1 - \exp ( -m / \tau_{t} )}{ ( m / \tau_{t} ) } ] + \beta_{2t} [ \frac{ 1 - \exp ( - m / \tau_{t} ) }{ (m / \tau_{t} ) } - \exp( - m / \tau_{t} ) ]$

ダイナミックなNelson-Siegelモデルにおいて，

β0は，カーブ全体の高低を規定する「水準ファクター（Level）」
β1は，カーブの角度を規定する「傾きファクター（Slope）」
β2は，カーブの湾曲の程度を規定する「曲率ファクター（Curvature）」

を表す。

　一方で今回の記事では，主成分分析を用いてイールドカーブの構造分解をおこなう*3。

主成分分析による推定では，イールドカーブの変動のほぼ全てが第一主成分から第三主成分の3つの共通成分によって説明できるとされている。

第一主成分は水準（Level），第二主成分は傾き（Slope），第三主成分は曲率（Curvature）を表していると解釈される。

データの取得

スポットレートのデータを使ってイールドカーブの構造を推定したかったが，データをうまくゲットできなかった（日本相互証券のHPとかいろいろ見たが）ので，財務省が公表している国債金利情報で代用*4。

取得するデータは，2001年1月から2017年7月までの月次データ（月末値）。年限は1年物から20年物まで。

プロットして動きを見てみる。 f:id:shiba_ryu0209:20170831104952j:plain

分析っぽいこと

国債金利データを直接，主成分分析をしてみる（定常性とか変化率云々の議論は知らない）。

読み込んだ金利データをxts型に変換し，princompを使って主成分分析。

#主成分分析
pca <- princomp((coredata(yield)), cor = T) #相関係数行列について主成分分析
summary(pca)

分析結果は，次のように。

Importance of components:
                          Comp.1     Comp.2     Comp.3     Comp.4       Comp.5       Comp.6       Comp.7       Comp.8       Comp.9      Comp.10      Comp.11
Standard deviation     3.2906484 0.96200235 0.46598691 0.12134893 0.0793525441 0.0574088587 0.0457253972 0.0381286418 0.0238068720 1.710747e-02 1.379945e-02
Proportion of Variance 0.9023639 0.07712071 0.01809532 0.00122713 0.0005247355 0.0002746481 0.0001742343 0.0001211494 0.0000472306 2.438879e-05 1.586874e-05
Cumulative Proportion  0.9023639 0.97948464 0.99757996 0.99880709 0.9993318240 0.9996064721 0.9997807064 0.9999018559 0.9999490865 9.999735e-01 9.999893e-01
                           Comp.12
Standard deviation     0.011308049
Proportion of Variance 0.000010656
Cumulative Proportion  1.000000000

第三主成分までで累積寄与率（Cumulative Proportion）が約99.76%となる。よって先述したように，イールドカーブの変動のほぼ全てが第三主成分までの3つの共通成分によって説明されるようだ。

国債の各年物に対するファクター・ローディングを見ると，次の通り。 f:id:shiba_ryu0209:20170831114836j:plain

また，分析したデータの期間におけるファクターの推移を見ると，次の通り。 f:id:shiba_ryu0209:20170831122037j:plain

基本的にイールドカーブの変動は第一主成分，すなわち水準ファクターによってもたらされたようだ。ただし，実際のイールドカーブの推移（フラットニング）に合わない印象。

-1を乗じると，現実に合うようないい具合になった↓。 f:id:shiba_ryu0209:20170831123824j:plain

補足・感想

証券アナリストジャーナルに主成分分析を用いてバーベルロング-ブレッドショート等の市場運用の観点から金利構造を分析している論文があったので，いつかそういう分析をしてみたい
米国債などでも分析してみたい
追記(2018.1.4)

どうやらカーブ取引など実務で使う際は利回りを事前に標準化してから，主成分分析にかけるようだ。この記事はそれを反映していないため…ちゃんとした推計とは言えないかもしれない…。

*1:たとえば，木島(1999)がある。

期間構造モデルと金利デリバティブ (シリーズ現代金融工学)

作者: 木島正明
出版社/メーカー: 朝倉書店
発売日: 1999/10/01
メディア: 単行本
購入: 2人クリック: 11回
この商品を含むブログ (1件) を見る

*2:http://www.math.ku.dk/~rolf/teaching/NelsonSiegel.pdf

*3:主成分分析は，永田・棟近(2001)を参照。

多変量解析法入門 (ライブラリ新数学大系)

作者: 永田靖,棟近雅彦
出版社/メーカー: サイエンス社
発売日: 2001/04
メディア: 単行本
購入: 2人クリック: 9回
この商品を含むブログ (1件) を見る

*4:http://www.mof.go.jp/jgbs/reference/interest_rate/index.htm

2017-08-28

RによるSVモデルの推計

概要

前回は日経平均株価のボラティリティ構造を，GARCH，EGARCHモデルで推計した

shiba-ryu0209.hatenablog.com

今回は日経平均株価のボラティリティ構造を，確率的分散変動（stochastic volatility；SV）モデルで推計してみた
Rパッケージstochvolの使い方をメモ

SVモデルについて

SVモデルの説明は，渡部(2000)とShepard(2005)*1を援用。

ボラティリティ変動モデル (シリーズ現代金融工学)

作者: 渡部敏明,木島正明
出版社/メーカー: 朝倉書店
発売日: 2000/06
メディア: 単行本
購入: 1人クリック: 17回
この商品を含むブログ (2件) を見る

株式のt期の収益率Rtは，t-1期の期待値Et-1(Rt)と予測不可能なショックεtで説明される。

$R_{t} = E_{t-1} ( R_{t} ) + \varepsilon_{t}$

予測不可能なショックεtは，変数σtと平均0，分散1のi.i.d系列にしたがう確率変数ztの積。

$\varepsilon_{t} = \sigma_{t} z_{t}, \ \sigma_{t} > 0, \ z_{t} ～ i.i.d., \ E [ z_{t} ] = 0, \ Var( z_{t} ) = 1$

σtの2乗を，ボラティリティと呼ぶ。

SVモデルは，ボラティリティの対数値の変動を線形ARMAによって定式化。通常，AR(1)モデル

$\ln ( \sigma_{t}^{2} ) = \omega + \phi \ln ( \sigma_{t-1}^{2} ) + \eta_{t}, \ \eta_{t}～i.i.d.N ( 0, \sigma_{ \eta }^{2} )$

が用いられる。Φが1に近ければ近いほど，ショックの持続性が高い。

データの取得

とりあえず，quantmodパッケージを使って日経平均株価の終値の日次データを取得。

（適当に）今回も前回同様にQQE導入以降で。ついでに対数差分をとって日次収益率に変換。

#パッケージの読込
library("quantmod")

#データの読込
nikkei <- getSymbols("^N225", src = "yahoo", from = as.Date("2013-04-04"),
                     to = as.Date("2017-08-18"), auto.assign = FALSE) 
nikkei.price <- nikkei$N225.Close

#対数差分をとって収益率に変換
nikkei.return <- diff(log(nikkei.price))[-1]
nikkei.return <- na.spline(nikkei.return)

得た日次収益率をプロットすると，次のように。 f:id:shiba_ryu0209:20170826195017j:plain

SVモデルの推計

Rパッケージstochvolを使って推計する。パッケージのマニュアル*2を見ると，

“Efficient algorithms for fully Bayesian estimation of stochastic volatility (SV) models via Markov chain Monte Carlo (MCMC) methods.”

と書いてある。

関数svsampleを使うと推定できる（らしい）ので，やってみる（事前分布は適当に設定）。

#SVモデルの推計
#サンプリング回数は10000回
res <- svsample(nikkei.return, draws = 10000, burnin = 2000,
                priormu = c(-10,1), priorphi = c(20,1.1), priorsigma = 0.1) 
summary(res, showlatent = FALSE)

推定結果は次の通り。

Summary of 10000 MCMC draws after a burn-in of 2000.

Prior distributions:
mu        ~ Normal(mean = -10, sd = 1)
(phi+1)/2 ~ Beta(a0 = 20, b0 = 1.1)
sigma^2   ~ 0.1 * Chisq(df = 1)

Posterior draws of parameters (thinning = 1):
            mean     sd      5%    50%    95%  ESS
mu        -8.971 0.1964 -9.2912 -8.968 -8.661 5809
phi        0.948 0.0169  0.9179  0.950  0.973  324
sigma      0.298 0.0448  0.2302  0.295  0.377  203
exp(mu/2)  0.011 0.0011  0.0096  0.011  0.013 5809
sigma^2    0.091 0.0276  0.0530  0.087  0.142  203

Φの推定値を見ると0.948であり，ボラティリティのショックには高い持続性があることがわかる。

パッケージ内のvolplot関数を使うと，推定ボラティリティの推移をプロットしてくれる。

volplot(res, col = c(2,1,2))

f:id:shiba_ryu0209:20170828153535p:plain

また，パッケージ内のparatraceplot関数を使うと，推定パラメータの推移をプロットしてくれる。

par(mfrow = c(3,1))
paratraceplot(res)

f:id:shiba_ryu0209:20170828154218p:plain

パッケージ内のparadensplot関数を使うと，事前分布と事後分布をプロットしてくれる。

par(mfrow = c(2,2))
paradensplot(res)

f:id:shiba_ryu0209:20170828154521p:plain

補足

推定ボラティリティの推移の図の横軸を推定期間の営業日にしたかったが，日付関連のプログラミングがよく分からなかったので今度挑戦したい（アホですみません）。
SVモデルの論文をいくつか読んで，事前分布の設定を勉強しないといけないなと実感
多変量モデルでも推計できるらしいので，いつか活用したい

*1:http://ora.ox.ac.uk/objects/uuid:cca137f7-e304-41a6-bee3-54de46b0a965

*2:http://cran.r-project.org/web/packages/stochvol/stochvol.pdf

2017-08-26

quantmodでデータを取得

概要

今更ながらquantmodパッケージの便利さに気付いたのでメモ
quantmodパッケージで簡単に金融データが取得可能

パッケージのインストール

Rパッケージquantmod内のgetSymbols関数を使うことで，金融関連のデータを簡単に取得可能。

quantmodのマニュアル*1を見ると， Yahoo!ファイナンス以外にもFREDなどから直接Rにデータを取り込むことが可能（らしい）。

まあ細かい部分は後で読むことにして，とりあえずインストール。

ついでにパッケージ内の関数getSymbolsを使って，適当に2017年以降の日経平均株価の日次データを取得。

#パッケージの読込
install.packages("quantmod")
library("quantmod")
#データの読込
nikkei <- getSymbols("^N225", src = "yahoo", from = as.Date("2017-01-01"),
                     to = as.Date("2017-08-26"), auto.assign = FALSE)

日経平均株価の毎営業日の，始値(Open)，高値(High)，安値(Low)，終値(Close)，出来高(Volume)，調整後終値(Adjusted)の6つが入ったxts型データを直接Rに取り込むことができる。

また，パッケージ内の関数chartSeriesを使うと，getSymbols関数で読み込んだデータでいい具合に図を描くことができる。

#図の描画
chartSeries(nikkei)

f:id:shiba_ryu0209:20170828123625p:plain

chartSeries関数の引数TAには，テクニカル分析のツールが入っている。試しにボリンジャー・バンドを入れてみる（4月からプロットしてみる）。

chartSeries(nikkei, subset = "2017-04::", TA = "addBBands()")

f:id:shiba_ryu0209:20170828130555p:plain

また，この他のプロット関数としてラインチャート（lineChart），バーチャート（barChart），ローソク足チャート（candleChart）がある。

たとえば，chandleChart関数を使ってみる。

candleChart(nikkei,theme = chartTheme("white.mono")) #背景を白色に
my.theme <- chartTheme(up.col = "blue", dn.col = "red", #上昇時に青色，
                       area = "white", bg.col = "white") #下落時に赤色のロウソク足チャート
candleChart(nikkei, theme = my.theme) #プロット

f:id:shiba_ryu0209:20170828125152p:plain

*1:https://cran.r-project.org/web/packages/quantmod/quantmod.pdf

2017-08-20

RによるGARCHモデル，EGARCHモデルの推計

概要

はてなブログを始めてみたので，試しに初投稿
日経225のボラティリティをGARCHモデル，EGARCHモデルで推計してみた
推計結果の解釈はともあれ，Rによる推計手順をメモ

日経225のデータの取得

　とりあえず，日銀によって現在行われている量的・質的金融緩和（QQE）以降の日経平均株価を取得。 QQEが導入されたのは2013年4月4日*1。

nikkei <- getSymbols("^N225", src = "yahoo", from = as.Date("2013-04-04"),
                     to = as.Date("2017-08-18"), auto.assign = FALSE)
nikkei.price <- nikkei$N225.Close

　日経平均株価の終値の推移は，次の図。 f:id:shiba_ryu0209:20170820144107j:plain

　終値の対数差分をとって，収益率の過程を計算。欠損値はスプライン補間。

nikkei.return <- diff(log(nikkei.price))[-1]
nikkei.return <- na.spline(nikkei.return)

f:id:shiba_ryu0209:20170820145225j:plain

GARCHモデルとEGARCHモデル

　モデルの説明は，渡部(2000)を援用。

ボラティリティ変動モデル (シリーズ現代金融工学)

作者: 渡部敏明,木島正明
出版社/メーカー: 朝倉書店
発売日: 2000/06/01
メディア: 単行本
購入: 1人クリック: 17回
この商品を含むブログ (2件) を見る

　 $GARCH(p,q)$ モデルは

$\varepsilon_{t} = R_{t} - E_{t-1} ( R_{t} ) \\ \varepsilon_{t} = \sigma_{t} z_{t}, \ \ \sigma_{t} > 0, \ \ z_{t}～i.i.d., \ \ E(z_{t}) = 0, \ \ Var(z_{t})= 1\\ \sigma_{t}^{2} = \omega + \sum_{i=1}^{p} \beta_{i} \sigma_{t-i}^{2} + \sum_{j=1}^{q} \alpha_{j} \varepsilon_{t-j}^{2}, \\ \omega > 0, \ \ \beta_{i}, \ \ \alpha_{j} > 0, \ \ i = 1,2,...,p; \ \ j = 1,2,...,q$

である。ただし， $R_t$ は株価の収益率， $\varepsilon_t$ は期待収益率の乖離（すなわち予測誤差）。 $z_t$ は正規分布である必要はない。

　 $EGARCH(p,q)$ モデルは

$\ln(\sigma_{t}^{2}) = \omega + \sum_{i=1}^{p} \beta \ln( \sigma_{t-i}^{2} ) + \sum_{j=1}^{q} \alpha_{j} [ \theta z_{t-i} + \gamma (|z_{t-j}| - E (|z_{t-j}|) ] \\ z_{t} = \varepsilon_t / \sigma_{t}$

である。 $z_t$ は予測誤差 $\varepsilon_t$ をボラティリティ $\sigma_t$ で除した変数。

　以下では，GARCHモデルとEGARCHモデルの次数 $p$ と $q$ はそれぞれ1とする*2。

　GARCHモデルとEGARCHモデルの違いは， $\gamma$ に表れる。現実のマーケットでは，株価の下落（負のショック）の方がマーケットに与えるインパクトが大きいといわれている（レバレッジ効果）。 GARCHモデルではショックが対称であるが，EGARCHモデルでは非対称にすることができる。

推計

　まずは，GARCHモデルの推計。次数 $(p,q)$ は $(1,1)$ 。推計はtseriesパッケージのgarch関数。

x <- garch(nikkei.return, order = c(1,1), trace = F)
summary(x)

推計結果は次のよう。

Call:
garch(x = nikkei.return, order = c(1, 1), trace = F)

Model:
GARCH(1,1)

Residuals:
     Min       1Q   Median       3Q      Max 
-5.20347 -0.46424  0.04545  0.59874  3.37857 

Coefficient(s):
    Estimate  Std. Error  t value Pr(>|t|)    
a0 7.942e-06   1.615e-06    4.917 8.79e-07 ***
a1 2.546e-01   3.130e-02    8.135 4.44e-16 ***
b1 7.560e-01   2.313e-02   32.678  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Diagnostic Tests:
    Jarque Bera Test

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


    Box-Ljung test

data:  Squared.Residuals
X-squared = 1.5052, df = 1, p-value = 0.2199

　係数はすべて1%水準で有意。β1は0.756であり，ボラティリティにはある程度の持続性があることがわかる。 Jarque Bara検定の結果をみると，p-valueが非常に小さいので残差が正規分布にしたがっていないと言える。 Box-Ljung検定の結果は，残差2乗に系列相関があるっぽいことを示している。

　次に，EGARCHモデルの推計。次数(p,q)は(1,1)。推計はrugarchパッケージ。

y <- ugarchspec(variance.model = list(model = "eGARCH", garchOrder = c(1,1)),
                mean.model = NULL,
                distribution.model = "norm") #今回はz_tに正規分布を仮定
y.fit = ugarchfit(y,nikkei.return)

　推計結果は次のよう。

*---------------------------------*
*          GARCH Model Fit        *
*---------------------------------*

Conditional Variance Dynamics   
-----------------------------------
GARCH Model : eGARCH(1,1)
Mean Model  : ARFIMA(1,0,1)
Distribution    : norm 

Optimal Parameters
------------------------------------
        Estimate  Std. Error   t value Pr(>|t|)
mu      0.000185    0.000339   0.54647  0.58474
ar1    -0.453453    0.028285 -16.03175  0.00000
ma1     0.412972    0.028990  14.24517  0.00000
omega  -0.339243    0.060716  -5.58738  0.00000
alpha1 -0.098662    0.017358  -5.68377  0.00000
beta1   0.959716    0.006908 138.91911  0.00000
gamma1  0.231997    0.029801   7.78492  0.00000

Robust Standard Errors:
        Estimate  Std. Error  t value Pr(>|t|)
mu      0.000185    0.000369   0.5016 0.615951
ar1    -0.453453    0.014511 -31.2485 0.000000
ma1     0.412972    0.014091  29.3083 0.000000
omega  -0.339243    0.070179  -4.8340 0.000001
alpha1 -0.098662    0.043259  -2.2807 0.022565
beta1   0.959716    0.007744 123.9298 0.000000
gamma1  0.231997    0.043629   5.3175 0.000000

LogLikelihood : 3174.769 

Information Criteria
------------------------------------
                    
Akaike       -5.8990
Bayes        -5.8666
Shibata      -5.8991
Hannan-Quinn -5.8867

Weighted Ljung-Box Test on Standardized Residuals
------------------------------------
                        statistic p-value
Lag[1]                     0.4827  0.4872
Lag[2*(p+q)+(p+q)-1][5]    1.6308  0.9946
Lag[4*(p+q)+(p+q)-1][9]    3.2827  0.8433
d.o.f=2
H0 : No serial correlation

Weighted Ljung-Box Test on Standardized Squared Residuals
------------------------------------
                        statistic p-value
Lag[1]                    0.01594  0.8995
Lag[2*(p+q)+(p+q)-1][5]   0.28856  0.9849
Lag[4*(p+q)+(p+q)-1][9]   0.74132  0.9945
d.o.f=2

Weighted ARCH LM Tests
------------------------------------
            Statistic Shape Scale P-Value
ARCH Lag[3] 5.569e-05 0.500 2.000  0.9940
ARCH Lag[5] 3.643e-01 1.440 1.667  0.9227
ARCH Lag[7] 3.928e-01 2.315 1.543  0.9870

Nyblom stability test
------------------------------------
Joint Statistic:  1.5595
Individual Statistics:              
mu     0.45356
ar1    0.10061
ma1    0.08621
omega  0.35353
alpha1 0.38347
beta1  0.36896
gamma1 0.09743

Asymptotic Critical Values (10% 5% 1%)
Joint Statistic:         1.69 1.9 2.35
Individual Statistic:    0.35 0.47 0.75

Sign Bias Test
------------------------------------
                   t-value   prob sig
Sign Bias          0.04305 0.9657    
Negative Sign Bias 0.07949 0.9367    
Positive Sign Bias 0.82657 0.4087    
Joint Effect       1.03200 0.7935    


Adjusted Pearson Goodness-of-Fit Test:
------------------------------------
  group statistic p-value(g-1)
1    20     59.07    5.420e-06
2    30     78.40    1.980e-06
3    40     90.69    5.337e-06
4    50    105.52    5.086e-06


Elapsed time : 0.3227599

　beta1は0.9597(1%有意)であり，ボラティリティには非常に高い持続性があることがわかる。また，gamma1は0.231997(1%有意)であり正である。このことから，QQE導入期に限れば日経平均株価のショック構造に非対称性がないことがわかる。

結果

　GARCHモデルとEGARCHモデルを使って，日経平均株価のボラティリティ構造をモデル化した。得られた結果は，

ボラティリティにはある程度以上の持続性があるらしい
ショック構造にはレバレッジ効果のような非対称性はないらしい

というもの。

補足

ショック部分に正規分布を仮定して推計したが，実際は正規分布より裾の厚い分布をしていることが指摘されている。なので，t分布や一般化誤差分布を仮定して推計した方がいいらしい（なおrugarchパッケージのugarchspec関数内の引数を変えることで，仮定する分布を変えることが可能）。
計算過程を確かめるためにパッケージの中身をちゃんと見ないとダメ。

*1:QQEについては，

http://www.boj.or.jp/mopo/outline/qqe.htm/

*2:次数は，天下り的に(1,1)に設定されることが多いらしい。詳しくは，Hansen and Lunde(2005),"Forecast Comparison of Volatility Models：Does Anything Beat A GARCH(1,1)？``http://onlinelibrary.wiley.com/doi/10.1002/jae.800/abstractや三井(2014)。

ARCH型モデルによる金融資産分析

作者: 三井秀俊
出版社/メーカー: 税務経理協会
発売日: 2014/12/11
メディア: 単行本
この商品を含むブログを見る

金融ネタとか備忘録とか

社会人6年目になりました（23年4月）。個人的に勉強したことをメモしていきます（2018年3月）。

日本国債金利の主成分分析

概要

Nelson-Siegelモデルについて

データの取得

分析っぽいこと

補足・感想

RによるSVモデルの推計

概要

SVモデルについて

データの取得

SVモデルの推計

補足

quantmodでデータを取得

概要

パッケージのインストール

RによるGARCHモデル，EGARCHモデルの推計

概要

日経225のデータの取得

GARCHモデルとEGARCHモデル

推計

結果

補足