【日次GDP】日本におけるGDPナウキャスティングの論文サーベイ

おはおんばんにちは。日次GDP予測やっていますが、なかなか当たらないですね。Giannoneたちの論文ではあんなに当たっていたのになぜなんでしょう。調べてみると、日本でもGDP等のナウキャスティングを実践する論文がいくつかあるようで、今回は先人たちの知恵にあやかろうということでサーベイを行いたいと思います。

① New Monthly Estimation Approach for Nowcasting GDP Growth: The Case of Japan（日銀ワーキングペーパー）

日本銀行の原さんと山根さんがによる研究で、約５００個の経済統計から月次GDPを推計する論文となっています。彼らのモデルの特徴は大きく２点です。１点目は経済統計を一気に主成分にかけるのではなく、鉱工業生産（製造業）と第3次産業活動指数（サービス業）だけは独立した説明変数としつつ、他の５００個近い経済統計は抽出した主成分説明変数として使用している点です。つまり、鉱工業生産（製造業）と第3次産業活動指数（サービス業）はGDPの大部分を説明するものとして、他の変数と独立して用いているということです。２点目は他の５００個の経済統計を需要項目、供給項目に分別し、また需要項目は消費、投資など細かいセクター分けをした後、それぞれをグループで主成分分析を行っている点です。また、GDP統計と同時点の説明変数を用いているのも特徴です。 具体的に見ていきましょう。まず、その他４７３系列の経済統計の中から鉱工業生産と第3次産業活動指数と相関のない部分を抽出します。以下の式をOLS推計します。

x_{t,m}^{i} = \alpha_0^{i}+\alpha_1^{i}dlog(IIP_{t,m})+\alpha_2^{i}dlog(ITA_{t,m})+\epsilon_{t,m}^{i}

ここで、x_{t,m}^{i}は四半期tの第mヶ月目の[tex;i]番目の経済変数になります（レベルのものはレベルのまま推計）。その結果、473個の\epsilon_{t,m}^{i}が推計できますが、これは鉱工業生産ならびに第3次産業活動指数とは相関のない部分となります（というかOLSの定義がそうなだけで実際相関がないとはいえない気がしますが）。これらの残差を主成分にかけるわけですが、先述の通りこの論文ではGDPの統計マニュアルに従い、これらを需要項目と供給項目に振り分け、それぞれに主成分をかけます。また、需要項目については、①消費、②投資、③国際貿易、④その他需要項目にセクター分けをし、供給と合わせて計５つのセクターでそれぞれ６つの主成分を抽出します。そして、以下の式をOLS推計します。

y_{t} = \beta_{0} + \beta_{1}dlog(IIP_{t}) + \beta_{2}dlog(ITA_{t}) + \beta_{3}p_{t}^{c} + \beta_{4}p_{t}^{i} + \beta_{5}p_{t}^{x} + \beta_{6}p_{t}^{o} + \beta_{7}p_{t}^{s} + \eta_{t}

ここで、y_{t}は季節調整済みの四半期GDP速報値、p_{t}^{z}は主成分であり、それぞれ消費、投資等上述したセクターを表しています。主成分は全部で３０個あるので、説明変数の選択としては6^{5}通りの組合せがあります。この論文では、AICが最も高くなる組合せを推計結果としています。ここで推計した係数を使って、月次GDPを推計していきます。推計式は四半期ベースなので、それを月次ベースに変換する必要があります。変換は以下の方程式を使用します。

y_{t} = \frac{1}{3}(z_{t,1}+z_{t,2}+z_{t,3})

ここで、

z_{t,1}=y_{t-1,2}+y_{t-1,3}+y_{t,1}

z_{t,1}=y_{t-1,3}+y_{t,1}+y_{t,2}

z_{t,1}=y_{t,1}+y_{t,2}+y_{t,3}

です。y_{t,m}は単月のGDP成長率、z_{t,m}は3カ月リンクしたGDP成長率の近似値を表しています。ここから、先ほどの推計式を月次版に修正したものを用いて、月次GDPを推計します。

\hat{y}_{t,m}=\frac{1}{3}\hat{\beta}_{0} + \hat{\beta}_{1}dlog(IIP_{t,m}) + \hat{\beta}_{2}dlog(ITA_{t,m}) + \hat{\beta}_{3}p_{t,m}^{c} + \hat{\beta}_{4}p_{t,m}^{i} + \hat{\beta}_{5}p_{t,m}^{x} + \hat{\beta}_{6}p_{t,m}^{o} + \hat{\beta}_{7}p_{t,m}^{s}

ハットがついているものは推計値を表しています。上式で計算した当該四半期の\hat{y}_{t,m}を平均することで、四半期のGDPを計算します。なお、月次鉱工業生産指数はラグがあるため、その予想値と日本版PMIにより計算した推計値を使用しています。また、第3次産業活動指数も同様に推計値を使用しています。 推計結果はまずまずで、もっともよいモデルでは自由度修正済み決定係数が0.921となっており、またリーマンショックや東日本大震災時の鋭い下落もとらえることが可能であるようです。また、エコノミストたちの予想値よりも的中率が高い結果も得られたようです。この論文では、二次速報値の推計も行っています。これまでの鉱工業生産と第3次産業活動指数に加え、法人企業統計の設備投資、投資項目、その他需要項目、供給項目の主成分を新たに説明変数に追加し、推計を試みています。その結果、法人企業統計の設備投資が速報値修正の大部分を担っていることがわかりました。また、こちらもエコノミストたちの予想より予測力が高い結果となりました。

② 景気判断における検索データの利用可能性（BOJ Reports & Research Papers）

日本銀行の白木さん、松村さん、松本さんによる論文です。こちらは経済統計ではなく、Google Trendで提供されているネット上の検索結果を用いて、観光庁「旅行業者取扱額」を説明しようとしています。このデータはGoogleの検索エンジンで任意の語句が検索された頻度を示す時系列データであり、週次データで2004年から公開されています（無料）。この論文では、Google Trendがセクター分けした「旅行カテゴリ」の語句（テーマパーク、ホテル、旅行代理店等）を月次データに変換し、季節調整をかけたものを使用しています。そして、それら21系列のデータを主成分分析で次元圧縮をしたものと景気ウォッチャー調査を説明変数とした以下の推計式で推計を行っています。

旅行取扱額（前期比）= \alpha + \beta 検索データ第１主成分１ + \gamma 景気ウォッチャー調査（前期差） + \epsilon

その結果、このモデルの予測力は景気ウォッチャー調査だけを説明変数として用いた推計式とあまり変わらないものの、リーマンショックや東日本大震災の鋭い落ち込みは良いことが示されています。

③国会議事録を用いた経済指標のナウキャスティング（DBSJ Japanese Journal）

早稲田大学の高杉さんと山名さんの論文です。この論文は上の二つと違い、テキストマイニング手法を用いて経済指標をナウキャストしています。具体的には、国会議事録の発言内容を一カ月ごとに分けたトレーニングデータを形態素解析にかけ、動詞、名詞、形容詞、副詞を抽出します。そして、それらをもとに一文単位でn-gramを作成します。なお、nをいくらにするかは被説明変数とする経済指標毎に最も予測力が高くなるよう探索的に調整をしています。最終的には、I種類のn-gramの出現回数を説明変数とするのですが、出現回数はトレーニングデータ毎に総和を取り、サンプル機関24カ月のパネルデータとしています（閾値を設定し、出現回数が極端に少ないものは除外）。次に、予測モデルの精度を上げるためにn-gramの数を絞ります。予測したい敬座指標とのピアソンの相関係数を計算し、閾値を超えるものに絞り込んだうえで、主成分分析を実行し、説明変数を絞り込みます。そうして絞り込んだ出現回数を説明変数、経済指標を被説明変数とした重回帰分析を行い、予測モデルを推定します。推計結果はまずまずといった感じで、AR(1)モデルによる推定よりも良いといった感じです。

いま行っている分析に近いのは１番目の論文ですので、ひとまずそちらの方向で分析を行っていきたいと思います。やはりもっと元データが必要ですね。また、個人消費や企業設備投資は二次速報で修正されることが多いので、検索データやテキストデータを用いて予測力の保管を行っていきたいと思います。

【備忘】VBAでネット上のファイルを自動でダウンロードする方法

備忘

おはこんばんにちは。 今日は備忘も備忘、VBAネタです。会社でVBAを使って、ファイルをダウンロードする方法について少し質問を受け、その回答に困ったので、ちょっとコードを書いてみたいと思います。

  • やりたいこと

「VBAを用いてDOM構造の中から欲しいファイルをタグ名を用いてダウンロードしたい」がやりたいことです。正直、Rを使ったりやSeleniumを併用したりするとすぐ取れるのですが、VBAはXpathにあまり対応してないらしく複雑なDOMには太刀打ちできない感じです。ご質問を受けた時は、この複雑なDOMから取得したい要素にどうやってたどり着けば良いのかとばかり考えていたのですが、ファイルをダウンロードするのが主たる目的ですのでそれならWindows APIを用いればそこまで難しくないかなと思いました。そのコードの実装を備忘として残しておきたいので、この記事を書いています。今回は実装の例として以前以下の記事で用いた財務省の金利情報のページから国債金利のcsvファイルをダウンロードするコードを書いてみたいと思います。

osashimix.hatenablog.com

  • VBAによる実装

まず、VBAでWEBスクレイピングを行うためには2つのライブラリファイルを参照可能にする必要があります。１つ目は「Windows HTML Object Library」で２つ目は「Microsoft Internet Controls」です。これはVBE画面のツール（T）タブから設定ができます（以下の画像参照）。この２つにチェックをつけて、OKを押下します。

f:id:osashimix:20181027112709p:plain

これで設定は完了です。財務省の金利情報のページはここです。このページの金利情報というリンクを押下すれば以下のようにファイルダウンロードの通知が現れます。

f:id:osashimix:20181027113534p:plain

今回はこのサイトのURLと要素を一意に絞るキーがわかっているとして、コードを書いていきます。今回の場合は以下の通り、jgbcm.csvがキーとなります。ファイルのダウンロードに関わる要素はリンク属性であるので、コードもキーを持つリンク属性のhrefを読み取り、それをURLDownloadToFile関数に渡すことでファイルをダウンロードする構成にします。

f:id:osashimix:20181027113901p:plain

では、コードを書いていきます。

' Winows APIのURLDownloadToFile関数を呼び込むためのおまじない
Declare Function URLDownloadToFile Lib "urlmon" Alias _
    "URLDownloadToFileA" (ByVal pCaller As Long, _
    ByVal szURL As String, _
    ByVal szFileName As String, _
    ByVal dwReserved As Long, _
    ByVal lpfnCB As Long) As Long

'ファイルダウンロードマクロを部品として作成
Sub fileDL(URL As String, _
           Key As String, _
           savepath As String)

Dim objIE As InternetExplorer
Dim fileName As String, link As String, saveName As String
Dim result As Long
Dim objLink As Object

' URLのホームページにアクセス
Set objIE = CreateObject("InternetExplorer.Application")
objIE.navigate URL

' ページの読み込み待ち
Do While objIE.Busy = True Or objIE.readyState < READYSTATE_COMPLETE 
 
    DoEvents
 
Loop

' DOMのリンク属性の中からKeyを含む要素を特定し、そのリンクを取得（一意に定まるKeyでないと特定できた最初のファイルをDL）
For Each objLink In objIE.document.Links

   If InStr(objLink.outerHTML, Key) > 0 Then

      link = objLink.href

      Exit For

   End If

 Next

If link = "" Then
    MsgBox "特定できませんでした"
    Exit Sub
  End If

' ファイル名とファイルを保存先を作成
fileName = Mid(link, InStrRev(linkValue, "/") + 1)
saveName = savepath & "\" & fileName

' ダウンロードを実行
result = URLDownloadToFile(0, link, saveName, 0, 0)

' 0の場合は成功、失敗の場合は-2146697210（resultをlongで定義する必要あり）
If result <> 0 Then
    MsgBox "ダウンロードできませんでした"
End If
End Sub

この部品をメインのマクロに組み込み使用すると、

Sub test()

    Call fileDL("https://www.mof.go.jp/jgbs/reference/interest_rate/", _
           "jgbcm.csv", _
           "C:\Users\aashi\Downloads")

End Sub

f:id:osashimix:20181027124602p:plain

無事、ダウンロードすることができました。やってみると、たいしたことない話ですね。ただ、一意に要素を絞るKeyがないとやはりVBAではWEBスクレイピングは難しそうです。こういうあたり、どうもVBAが好きになれないところです。会社でもPythonやRが使えたらどれだけ作業効率が上がることか。まあとりあえず、当面はVBAでどうにかこうにかやるしかなさそうです。今週はできればもう一本投稿したいと思っているので、そちらもご笑覧ください。 では。

【仕事関連】イールドカーブのモデル化について

その他

お久しぶりです。実はPCがウイルスにかかってしまい、新しいPCを買うまでにお金をためていたので、なかなかブログを更新できていませんでした。今回は、今やっている投資顧問会社でのレポーティングに仕事に関連する投稿です。債券運用のパフォーマンス要因分解を行う際、①イールドカーブ効果、②銘柄選択効果、③為替効果で超過リターンを分解することがあります。そして、このイールドカーブ効果は理論イールドカーブを使用して算出するらしいのですが、イールドカーブってどうやってモデル化するのだろうと思ったので、今回まとめてみたいと思います。ぱっと浮かぶのはやはりNelson-Siegelモデルです。マクロ経済予測なんかでも使用されているモデルです。まずこの復習から始めたいと思います。

Nelson-Siegelモデル

NelsonとSiegelが1987にJournal of Businessにパブリッシュした論文です。http://www.jstor.org/stable/2352957

このモデルの特徴は少ないパラメータでイールドカーブの一般的な特徴を記述できるところです。このモデルでは、満期mヵ月の無リスクフォワードレートf(m)を

{ \displaystyle f(m) = L + Se^{-m\lambda} + Cm\lambda e^{-m\lambda}} \tag{1}

と表しています。ここで、L、S、C、\lambdaはパラメータであり、それぞれイールドの水準、傾き、曲率、スケールを決定します。y(m)を残存mヵ月の無リスク割引債のスポットレートだとすると、y(m)は以下のようなフォワードレートf(m)の平均金利として表せます。

\displaystyle y(m) = \frac{1}{m} \int_0^m f(x) dx \tag{2}

\displaystyle = \frac{1}{m} \int_0^m (L + Se^{-x\lambda} + Cx\lambda e^{-x\lambda}) dx \tag{3}

\displaystyle = \frac{1}{m} (\int_0^m L dx + \int_0^m Se^{-x\lambda} dx + \int_0^m Cx\lambda e^{-x\lambda} dx) \tag{4}

各積分はそれぞれ以下のように計算できます。

\displaystyle \int_0^m L dx = Lm

\displaystyle \int_0^m Se^{-x\lambda} dx = S[-\frac{e^{-x\lambda}}{\lambda}]_0^m = S(\frac{1}{\lambda} - \frac{e^{-m\lambda}}{\lambda})

\displaystyle \int_0^m Cx\lambda e^{-x\lambda} dx = C\int_0^m (x\lambda e^{-x\lambda} - e^{-x\lambda} + e^{-x\lambda}) dx = C[-\frac{e^{-x\lambda}}{\lambda} -xe^{-x\lambda}]_0^m = C(\frac{1}{\lambda} - \frac{e^{-m\lambda}}{\lambda} - me^{-m\lambda})

よって、(4)は、

\displaystyle y(m) = \frac{1}{m} (Lm + S(\frac{1-e^{-m\lambda}}{\lambda}) + C(\frac{1-e^{-m\lambda}}{\lambda} - me^{-m\lambda})) \tag{5}

と表すことができます。先ほど少しパラメータの意味についてご説明しましたが、この形を見ればよりその意味が分かるのではないかと思います。(5)の第一項は定数項なので、全ての残存期間に対して共通の金利水準を表しています（つまり、金利シフトに関わるパラメータ）。第二項は残存期間に対して以下のように単調減少関数となっており、Sはイールドカーブの勾配を表していることがわかります。つまり、おおざっぱに言えば、 S>0の時は逆イールド、逆の時は順イールドを表現することができるということです。

f:id:osashimix:20181007121721p:plain

次に、第三項ですがこちらは残存期間に対して上に凸の関数であり、Cは曲率を表していることがわかります。

f:id:osashimix:20181007122945p:plain

これらの項を組み合わせると、以下のように見慣れた順イールドカーブを表現することができ、また先述したようにパラメータの値によっては逆イールドも表現でき、また勾配と曲率の組合せによりフラットニングやスティープニングも表現できます。*1

f:id:osashimix:20181007021925p:plain

あとは(5)を実際のゼロクーポンを当ててフィッティングすれば良いことになります。ここで、L、S、Cに関しては線形なので通常のOLSで推計すればよいことになりますが、\lambdaは非線形関数となっており、optimizerによる収束計算が必要になります。

Rによる推定の実装

推定を行うにあたって必要なものはゼロクーポンイールドですが、今回は財務省のデータ（ 国債金利情報 : 財務省）で代用したいと思います。*2プロットするとこんな感じ（2000年4月~直近）。

f:id:osashimix:20181008012745p:plain

金利は低下傾向であり、また2010年の前と後では短期ゾーンの形状が変わっていることもわかると思います。つまり、これは先述したパラメータが一定ではないことを意味しており、各時期のイールドカーブにフィッティングするためには時変パラメータの推定が必要であることを示唆しています。当たり前といえば当たり前ですが、しかしパラメータが推定期間に対して安定性がないのではモデルとしてはいまいちなのではと思ってしまいます。ちなみにイールドカーブのコードは以下です。

library(dplyr)
library(tidyr)
library(ggplot2)
library(scales)
library(xts)
library(plotly)

source.jgb <- "http://www.mof.go.jp/jgbs/reference/interest_rate/data/jgbcm_all.csv"
jgb_all <- read.csv(source.jgb,skip=1)
colnames(jgb_all) <- c("Date","1","2","3","4","5","6","7","8","9","10",
                       "15","20","25","30","40")
jgb_all$Date <- as.character(jgb_all$Date)

# a function to transform Japanese calendar to Western calendar 
ToChristianEra <- function(x)
{
  era  <- substr(x, 1, 1)
  year <- as.numeric(substr(x, 2, nchar(x)))
  if(era == "H"){
    year <- year + 1988
  }else if(era == "S"){
    year <- year + 1925
  }
  as.character(year)
}

# calendar transformation (Japanese to Western)
jgb.day <- strsplit(jgb_all[,1],"\\.")
warn.old <- getOption("warn")
options(warn = -1)
jgb.day <- lapply(jgb.day, function(x)c(ToChristianEra(x[1]), x[2:length(x)]))
jgb_all[, 1] <- as.Date(sapply(jgb.day, function(x)Reduce(function(y,z)paste(y,z, sep="-"),x)))
jgb_all[, -1] <- apply(jgb_all[, -1], 2, as.numeric)
options(warn = warn.old)
rm(jgb.day)

# create a xts object
jgb.xts <- as.xts(read.zoo(jgb_all))

# 3D plot
jgb.xts["2000::"] %>%
  data.matrix() %>% 
  t() %>%
  plot_ly(
    x=as.Date(index(jgb.xts["2000::"])),
    y=c(1,2,3,4,5,6,7,8,9,10,15,20,25,30,40),
    z=.,
    type="surface"
  ) %>%
  plotly::layout(
    scene=list(
      xaxis=list(title="date"),
      yaxis=list(title="term"),
      zaxis=list(title="yield")
    )
  )

では、固定パラメータの推計に移ります。Nelson and Siegel(1987)では\lambdaの値を先に決め、その後でほかのパラメータをOLS推計することを繰り返し、決定係数が最大になる点を推定値としています。とりあえず、こちらもこの推計方法に従いたいと思います。*3まず推計用のデータセットを作成します。

# create dataset
dataset <- jgb_all %>%
  tidyr::gather(key = "maturity", value = "irate",
                "1","2","3","4","5","6","7","8","9","10","15","20","25","30","40")

> head(dataset)
        Date maturity  irate    
1 1974-09-24        1 10.327 
2 1974-09-25        1 10.333
3 1974-09-26        1 10.340
4 1974-09-27        1 10.347
5 1974-09-28        1 10.354 
6 1974-09-30        1 10.368

次に、決定係数を最大化するために目的関数を定義します。この目的関数は\lambdaを入力値とし、その値に応じて(5)式の第二項、第三項を計算します。そして、それらを説明変数、財務省から取得したイールドデータを非説明変数とする重回帰を実施し、その決定係数を返します。

# define objective function
obj.func <- function(lambda,dataset,time_end,time_start){
  f1 <- function(m,lambda){
    ((1-exp(-m*lambda))/(m*lambda))
    }
  f2 <- function(m,lambda){
    ((1-exp(-m*lambda))/(m*lambda)-exp(-m*lambda))
    }
  dataset <- dataset %>%
    mutate(f1 = f1(as.numeric(dataset$maturity)*12,lambda), f2 = f2(as.numeric(dataset$maturity)*12,lambda))
  results <- lm(0.01*irate ~ f1 + f2,data = dataset,subset = (Date>time_start & time_end>Date))
  return(summary(results)$r.squared)
}

Rのoptimaze関数で最適化を行います。サンプル自体は1974年からありますが、国債流通市場が機能し始めたのは1980年代後半の金利自由化以降であり、また残存期間10年超の国債が安定的に取引されていたのは1990年以降であるので、1992年以降のサンプルを使用して推計を行います（それまではNAがたくさん）。\lambdaはNelson and Siegel(1987)と同じく1/200~1/10の間で探索を行います。

optimize(obj.func,interval = c(1/10,1/200),dataset=dataset,time_start="1992-01-01",time_end="2018-09-28",maximum = TRUE)
$`maximum`
[1] 0.06396279

$objective
[1] 0.1468183

決定係数を最大にする\lambdaは0.064であり、その時の決定係数は0.147という結果になりました。いや、決定係数低すぎるなぁ。やはり、先ほどお示ししたように2010年（もしかするとリーマン）の前後でイールドカーブの形状が大きく変わっており、また全ての残存期間で金利が低下傾向ということが原因であると思います。やはり、固定パラメータでは全サンプルに対する当てはまりが悪い。安直な手ですが、サンプル期間を2010年で分断し、推計を二回行ってみます。もちろん、決定係数にこだわりすぎることは危険ではあるのですが、これではモデルを使用している意味がないほど当てはまりが悪いので。

# 1992-01-01~2009-12-31で推計
optimize(obj.func,interval = c(1/10,1/200),dataset=dataset,time_end="2010-01-01",time_start="1992-01-01",maximum = TRUE)
$`maximum`
[1] 0.06885455

$objective
[1] 0.1891524

# 2010-01-01~2018-09-28で推計
optimize(obj.func,interval = c(1/10,1/200),dataset=dataset,time_end="2018-09-28",time_start="2010-01-01",maximum = TRUE)
$`maximum`
[1] 0.02323009

$objective
[1] 0.6573314

どうやら2010年以降のサンプルでは当てはまり良いようです。この\lambdaの値を用いて、OLS推計を行います。

lambda = 0.02323009
dataset <- dataset %>%
  mutate(f1 = f1(as.numeric(dataset$maturity)*12,lambda), f2 = f2(as.numeric(dataset$maturity)*12,lambda))
> head(dataset)
        Date maturity  irate        f1        f2
1 1974-09-24        1 10.327 0.8727162 0.1159956
2 1974-09-25        1 10.333 0.8727162 0.1159956
3 1974-09-26        1 10.340 0.8727162 0.1159956
4 1974-09-27        1 10.347 0.8727162 0.1159956
5 1974-09-28        1 10.354 0.8727162 0.1159956
6 1974-09-30        1 10.368 0.8727162 0.1159956

results <- lm(0.01*irate ~ f1 + f2,data = dataset,subset = (Date>"2010-01-01" & "2018-09-28">Date))
> summary(results)

Call:
lm(formula = 0.01 * irate ~ f1 + f2, data = dataset, subset = (Date > 
    "2010-01-01" & "2018-09-28" > Date))

Residuals:
       Min         1Q     Median         3Q        Max 
-0.0156039 -0.0024211  0.0005173  0.0025148  0.0098866 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.0207849  0.0000752  276.39   <2e-16 ***
f1          -0.0193824  0.0001003 -193.21   <2e-16 ***
f2          -0.0309228  0.0003217  -96.12   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.004046 on 32172 degrees of freedom
Multiple R-squared:  0.6573, Adjusted R-squared:  0.6573 
F-statistic: 3.086e+04 on 2 and 32172 DF,  p-value: < 2.2e-16

推計したイールドカーブをプロットしたのが以下です。 f:id:osashimix:20181008184635p:plain んー、なんか違う気がするなぁ。まあとりあえず、イールドカーブの推計はできました。

パフォーマンス要因分解とイールドカーブモデル

そもそもイールドカーブモデルはパフォーマンス要因分解のどこで使用するのでしょうか？以下のような参考文献を見つけました。p22のAppendix2にその説明があります。

http://summit.sfu.ca/system/files/iritems1/762/FRM%202010%20Melnikov,%20A.%20Simic,%20S..pdf

FIXED INCOME ATTRIBUTION: ANALYZING SOURCES OF RETURN Appendix2: DMT (at t-1) is a treasury rate on a 2.14 year treasury bill at the beginning of the attribution period. DMT (at t) is a treasury rate on a 2.14 year treasury bill at the end of the attribution period. It is unlikely that we are going to find 2.14 year treasury bill trading in the market at any given point in time. As such, we will be required to interpolate it’s yield from a standard treasury yield curve. There are several choices available for interpolation, with the simplest one being linear interpolation. Models that are more complex may apply quadratic, cubic interpolation, or Nelson-Siegel (1987) approach. As long as interpolation approach is consistent for both benchmark and portfolio, the bias is kept to minimum.

トータルリターンをシフト効果とツイスト効果に分解する際に、Duration-matched Tresury rateを計算する必要があるらしく、その際にイールドカーブモデルが必要になるとのことです。保有している債券のデュレーションはパフォーマンス計測期間にきっちり整数値をとるとは限りませんが、市場には例えば2.14のデュレーションを持つ債券は常にありそうにはないので当然です。ただ、Nelson Siegelモデルを使用するほかにも、線形補間や2次、3次補間を行うこともあるようです。こうして考えてみると、このモデルをパフォーマンス要因分解が目的で使用するのであれば、推計時点でパラメータの安定性がないことにさほど問題はないように思われます。時系列分析に使用するのが問題のような気がしてきました（いくつかペーパーがありますが）。基本的に補間のために作られたモデルと理解した方が良さそうです。Appendix3では各要因分解の計算方法が掲載されているので、それを見て今回の記事は終わりにしたいと思います。トータルリターンは以下のように、①Income、②Shift、③Twist、④Spread、⑤Selectionに分解することができます。

Income = \frac{Coupon}{Beginning Price}

Shift = -Duration\times\Delta KRD *4

Twist = -Duration\times (\Delta DMT - \Delta KRD)

Spread = -Duration\times (\frac{BM Total Return - BM Income - BM Shift - BM wist}{BM Duration})

Selection = Total Return - Income - Shift - Twist - Spread

*1:Nelson and Siegel(1987)では、L、S、Cをそれぞれ長期、短期、中期的な影響を与えるパラメータと解釈しています。その理由として、Lが残存期間に対して不変でゼロにはならないこと、Sにかかる項が他の項と比較して急速、単調に０へ減退すること、Cにかかる項が０から始まり、（極限では）０に収束することを挙げています。

*2:1年～40年ま でを1年毎に区切り、3次ス プライン関数を用いて補間することにより、コン スタントマチュリティーベースの金利（年限が常 に一定である金利）を算出しています。

*3:Nelson and Siegel(1987)では各時点で得られたイールドカーブデータに対して個別に推計を行っていますが、ここでは全サンプル一気に推計をしています。つまり、Nelsonらは\lambdaを時変パラメータのように扱っていますが、こちらでは全てのサンプル期間で一定としています。

*4:\Delta KRDは5年または10年のkey rate durationの対前期変化

【備忘】Black-Scholes-Mertonの公式の証明

備忘

今回は完全に自己満足な記事となっておりますので、その点ご留意ください。

おはこんばんにちは。運用部ではないにしろ投資顧問会社で働いているとデルタやセータなどのグリークスに出会うことがあります。そんな時、そもそもブラックショールズ方程式ってどうやって証明するんだっけと思うことが良くありました。そこで、この記事ではジョン・ハルのフィナンシャルエンジニアリング（第7版）の第13章付録にあるBlack-Scholes-Mertonの公式の証明を復習したいと思います。

フィナンシャルエンジニアリング―デリバティブ取引とリスク管理の総体系

フィナンシャルエンジニアリング―デリバティブ取引とリスク管理の総体系

補題の証明

まず証明に必要な補題の証明を行いたいと思います。

（ジョン・ハル「フィナンシャルエンジニアリング（第7版）」P453より）

Vが対数正規分布に従い、lnVの標準偏差がwであるとき、 E[max(V-K,0)]=E(V)N(d_{1})-KN(d_{2})\tag{13A.1} が成り立つ。ここで、

d_{1} = \frac{ln[E(V)/K]+w^{2}/2}{w}

d_{2} = \frac{ln[E(V)/K]-w^{2}/2}{w}

であり、Eは期待値を表す。（引用終了）

なお、Kは定数です。

証明

g(V)をVの確率密度関数とします。このとき、期待値の定義から \displaystyle E[max(V-K,0)]=\int_K^{\infty} (V-K)g(V)dV\tag{13A.2} となります。なお、右辺の積分区間がKから\inftyとなっているのは、K\leqq0でmax(V-K,0)は0であるからです（書いてもいいけど省略している）。Vは対数正規分布に従うので、その期待値mは、 m=ln[E(V)]-\frac{w^{2}}{2}\tag{13A.3} となります。ここで、Qを、 Q=\frac{lnV-m}{w}\tag{13A.4} と定義すると、この変数は期待値0、標準偏差1.0の標準正規分布に従います（Vは対数正規分布に従うのでlnVは正規分布に従う）。要は対数変換したのち、正規化したということです。なので、密度関数h(Q)は、 h(Q)=\frac{1}{\sqrt{2\pi}}\exp(-Q^{2}/2) となります。（13A.4）を用いて、（13A.2）式の右辺の積分でVをQに変数変換すると、以下のようになります。まず、（13A.4）をlnVについて解くと、

lnV=Qw+m

となり、ここから

V=e^{Qw+m}

となります。また、（13A.4）をVについて微分し、dVについて解くと、

dV=wVdQ

が得られます。Vは対数正規分布に従うので、その確率密度関数g(V)は、

g(V)=\frac{1}{\sqrt{2\pi}w V}exp\{-\frac{(lnV-m)^{2}}{2w^{2}}\}

なので、これらを（13A.2）に代入し、

\displaystyle E[max(V-K,0)]=\int_{(lnK-m)/w}^{\infty} (e^{Qw+m}-K)\frac{1}{\sqrt{2\pi}w V}exp\{-\frac{(Qw+m-m)^{2}}{2w^{2}}\}wVdQ

整理すると、

\displaystyle E[max(V-K,0)]=\int_{(lnK-m)/w}^{\infty} (e^{Qw+m}-K)h(Q)dQ

が得られます。上式を書き換えると、 \displaystyle E[max(V-K,0)]=\int_{(lnK-m)/w}^{\infty} e^{Qw+m}h(Q)dQ-K\int_{(lnK-m)/w}^{\infty}h(Q)dQ\tag{13A.5} が得られます。ここまではmax(V-K,0)の期待値を変数変換しただけです。ここで、

e^{Qw+m}h(Q)=\frac{1}{\sqrt{2\pi}}e^{(-Q^{2}+2Qw+2w)/2} \\ =\frac{1}{\sqrt{2\pi}}e^{[-(Q-w)^{2}+2m+w^{2}]/2}\\=\frac{e^{m+w^{2}/2}}{\sqrt{2\pi}}e^{[-(Q-w)^{2}]/2}\\=e^{[-(Q-w)^{2}]/2}h(Q-w)

となるので、(13A.5)式は

\displaystyle E[max(V-K,0)]=e^{m+w^{2}/2}\int_{(lnK-m)/w}^{\infty} h(Q-w)dQ \\ -K\int_{(lnK-m)/w}^{\infty}h(Q)dQ\tag{13A.6}

と書き換えられます。N(x)を標準正規分布の累積密度関数とすると、(13A.6)式の最初の積分はQが(lnK-m)/w-w以上になる確率なので、

1-N[(lnK-m)/w-w]

つまり、正規分布の対称性から

N[(-lnK+m)/w+w]

と表すことができます。(13A.3)式をmに代入すると、

N(\frac{ln[E(V)/K]+w^{2}/2}{w})=N(d_{1})

になります。同様に(13A.6)式の２番目の積分はN(d_{2})となる。よって、(13A.6)式は

E[max(V-K,0)]=e^{m+w^{2}/2}N(d{1})-KN(d_{2})

となり、(13A.3)式を代入すると証明完了。

Black-Scholes-Mertonの結果

ここまでこれば8割がた証明は終わっています。時点Tに満期を迎える配当のない株式に対するコールオプションを考えます。行使価格をK、無リスク金利をr、現在の株価をS_{0}、ボラティリティを\sigmaとします。コール価格cはこのコールオプションで得られる便益を無リスク金利で割り引いた値となります。

c=e^{-rT}E[max(S_{T}-K,0)]

ここで、S_{T}は時点Tにおける株価で対数正規分布に従います。Eはリスク中立世界における期待値を表しています。さらに、E(S_{T})=S_{0}e^{rT} *1、lnS_{T}の標準偏差は\sigma \sqrt{T}となります。先ほどの補題より（13A.7）式は、

c=e^{-rT}[S_{0}e^{rT}N(d_{1})-KN(d_{2})]\\=S_{0}N(d_{1})-Ke^{-rT}N(d_{2})

ここで、

\displaystyle d_{1}=\frac{ln[E(S_{T})/K]+\sigma^{2}T/2}{\sigma\sqrt{T}}=\frac{ln(S_{0}/K)+(r+\sigma^{2}/2)T}{\sigma\sqrt{T}}

\displaystyle d_{2}=\frac{ln[E(S_{T}/K)]-\sigma^{2}T/2}{\sigma\sqrt{T}}=\frac{ln(S_{0}/K)+(r-\sigma^{2}/2)T}{\sigma\sqrt{T}}

これで、Black-Scholes-Mertonの公式を導出できました。

*1:リスク中立世界においては、すべての人はリスクに関して無差別であるから、投資家はリスクに対する見返りを要求しない。よって、全ての証券の期待収益率は無リスク金利となる。

【競馬】【再実行】rvestでyahoo競馬にある過去のレース結果をスクレイピングしてみた。

今回で3回目になりますが、またまたrvestで過去のレース結果を落としてみたいと思います。過去の記事を見てないという人は先にそちらをご覧になられることをお勧めします。

osashimix.hatenablog.com

osashimix.hatenablog.com

今回データを取り直そうと思ったのは、競馬の分析をした際により多くの項目を説明変数に加えて、分析をしたいと思ったからです。なので、今回は前回のRスクリプトに追記を行う形でプログラムを作成しました（前回も前々回のスクリプトに付け足しただけですが）。新たに追加したデータ項目は以下の14個です。

  1. 芝かダートか
  2. 右回りか左回りか
  3. レースコンディション（良や稍重など）
  4. 天候
  5. 馬の毛色（栗毛、鹿毛など）
  6. 馬主
  7. 生産者
  8. 産地
  9. 生年月日
  10. 父馬
  11. 母馬
  12. そのレースまでの獲得賞金（2003年から入手可能）
  13. ジョッキーの体重
  14. ジョッキーの体重の増減

実はまだデータ収集は終わっていなくて、Rのプログラムがずっと実行中になっています（3日くらい回しています）。しかし、プログラム自体はきっちり回っているのでスクリプトの紹介をしていこうと思います。もしかしたら追記で結果を書くかもしれません。

スクリプトの中身

まずはパッケージの呼び出しです。

# rvestによる競馬データのwebスクレイピング

#install.packages("rvest")
#if (!require("pacman")) install.packages("pacman")
#install.packages("beepr")
#install.packages("RSQLite")
pacman::p_load(qdapRegex)
library(rvest)
library(stringr)
library(dplyr)
library(beepr)
library(RSQLite)

かなりwarnningが出るのでそれを禁止し、SQLiteに接続しています

# warnning禁止
options(warn=-1)

# SQLiteへの接続
con = dbConnect(SQLite(), "horse_data.db", synchronous="off")

1994年からしかオッズが取れないので、1994年から直近までのデータを取得します。yahoo競馬では月ごとにレースがまとめられているので、それを変数として使用しながらデータをとっていきます。基本的には、該当年、該当月のレース結果一覧へアクセスし、そのページ上の各日の個々の競馬場ごとのタイムテーブルへのリンクを取得します。個々の競馬場でレースはだいたい12ほどあるので、そのリンクを取得し、各レースのレース結果ページにアクセスします。そして、レース結果を取得していきます。まず、各日の個々の競馬場ごとのタイムテーブルへのリンクの取得方法です。

for(year in 1994:2018){
  start.time <- Sys.time() # 計算時間を図る
  # yahoo競馬のレース結果一覧ページの取得
  for (k in 1:12){ # kは月を表す
    
    tryCatch(
      {
        keiba.yahoo <- read_html(str_c("https://keiba.yahoo.co.jp/schedule/list/", year,"/?month=",k)) # 該当年、該当月のレース結果一覧にアクセス
        Sys.sleep(2)
        race_lists <- keiba.yahoo %>%
          html_nodes("a") %>% 
          html_attr("href") # 全urlを取得
        
        # 競馬場ごとの各日のレースリストを取得
        race_lists <- race_lists[str_detect(race_lists, pattern="race/list/\\d+/")==1] # 「result」が含まれるurlを抽出
      }
      , error = function(e){signal <- 1}
    )

ここでは、取得したリンクのurlにresultという文字が含まれているものだけを抽出しています。要はそれが各競馬場のレーステーブルへのリンクとなります。ここからは取得した競馬場のレーステーブルのリンクを用いて、そのページにアクセスし、全12レースそれぞれのレース結果が掲載されているページのリンクを取得していきます。

    for (j in 1:length(race_lists)){ # jは当該年月にあったレーステーブルへのリンクを表す
      
      tryCatch(
        {
          race_list <- read_html(str_c("https://keiba.yahoo.co.jp",race_lists[j]))
          race_url <- race_list %>% html_nodes("a") %>% html_attr("href") # 全urlを取得
          
          # レース結果のurlを取得
          race_url <- race_url[str_detect(race_url, pattern="result")==1] # 「result」が含まれるurlを抽出
        }
        , error = function(e){signal <- 1}
      )

各レース結果へのリンクが取得できたので、ここからはいよいよレース結果の取得とその整形パートに入ります。かなり長ったらしく複雑なコードになってしまいました。レース結果は以下のようなテーブル属性に格納されているので、まずそれを単純に引っ張ってきます。 f:id:osashimix:20180814213644p:plain

      for (i in 1:length(race_url)){ # iは当該年月当該競馬場で開催されたレースを表す
        
        print(str_c("現在、", year, "年", k, "月",j, "グループ、", i,"番目のレースの保存中です"))
        
        tryCatch(
          {
            race1 <-  read_html(str_c("https://keiba.yahoo.co.jp",race_url[i])) # レース結果のurlを取得
            signal <- 0
            Sys.sleep(2)
          }
          , error = function(e){signal <- 1}
        )
        
        # レースが中止orこれまでの過程でエラーでなければ処理を実行
        if (identical(race1 %>%
                      html_nodes(xpath = "//div[@class = 'resultAtt mgnBL fntSS']") %>%
                      html_text(),character(0)) == TRUE && signal == 0){
          
          # レース結果をスクレイピング
          race_result <- race1 %>% html_nodes(xpath = "//table[@id = 'raceScore']") %>%
            html_table()
          race_result <- do.call("data.frame",race_result) # リストをデータフレームに変更
          
          colnames(race_result) <- c("order","frame_number","horse_number","horse_name/age","time/margin","passing_rank/last_3F","jockey/weight","popularity/odds","trainer") #　列名変更

tableをただ取得しただけでは以下のように\nが入っていたり、一つのセルに複数の情報が入っていたりと分析には使えないデータとなっています。

> head(race_result)
      order frame_number horse_number                         horse_name/age     time/margin
    1     1            4            5     ヤマカツエース\n    \n牡4/492(+6)/          2.01.2
    2     2            5            7   マイネルフロスト\n    \n牡5/488(+4)/   2.01.33/4馬身
    3     3            6           10         フルーキー\n    \n牡6/490(+6)/   2.01.43/4馬身
    4     4            8           14 ライズトゥフェイム\n    \n牡6/488(+4)/ 2.01.61 1/4馬身
    5     5            3            3     ステラウインド\n    \n牡7/490(+4)/    2.01.6アタマ
    6     6            4            6 ブライトエンブレム\n    \n牡4/492(+4)/    2.01.6アタマ
      passing_rank/last_3F         jockey/weight popularity/odds   trainer
    1      05-05-04-0233.0   池添 謙一\n    56.0      3    (4.9) 池添 兼雄
    2      01-01-01-0134.4   松岡 正海\n    57.0      5    (6.9)   高木 登
    3      08-09-08-0832.7  M.デムーロ\n    57.5      1    (3.4) 角居 勝彦
    4      10-11-10-1032.6 石川 裕紀人\n    56.0     6    (12.4) 加藤 征弘
    5      03-03-03-0233.6   蛯名 正義\n    56.0     7    (17.2) 尾関 知人
    6      06-06-06-0633.2  C.ルメール\n    56.0      2    (4.8) 小島 茂之

なので、これを成型する必要が出てきます。

          # 通過順位と上り3Fのタイム
          race_result <- dplyr::mutate(race_result,passing_rank=as.character(str_extract_all(race_result$`passing_rank/last_3F`,"(\\d{2}-\\d{2}-\\d{2}-\\d{2})|(\\d{2}-\\d{2}-\\d{2})|(\\d{2}-\\d{2})")))
          race_result <- dplyr::mutate(race_result,last_3F=as.character(str_extract_all(race_result$`passing_rank/last_3F`,"\\d{2}\\.\\d")))
          race_result <- race_result[-6]
          
          # タイムと着差
          race_result <- dplyr::mutate(race_result,time=as.character(str_extract_all(race_result$`time/margin`,"\\d\\.\\d{2}\\.\\d|\\d{2}\\.\\d")))
          race_result <- dplyr::mutate(race_result,margin=as.character(str_extract_all(race_result$`time/margin`,"./.馬身|.馬身|.[:space:]./.馬身|[ア-ン-]+")))
          race_result$margin[race_result$order==1] <- "トップ"
          race_result$margin[race_result$margin=="character(0)"] <- "大差"
          race_result$margin[race_result$order==0] <- NA
          race_result <- race_result[-5]
          
          # 馬名、馬齢、馬体重
          race_result <- dplyr::mutate(race_result,horse_name=as.character(str_extract_all(race_result$`horse_name/age`,"[ァ-ヴー・]+")))
          race_result <- dplyr::mutate(race_result,horse_age=as.character(str_extract_all(race_result$`horse_name/age`,"牡\\d+|牝\\d+|せん\\d+")))
          race_result$horse_sex <- str_extract(race_result$horse_age, pattern = "牡|牝|せん")
          race_result$horse_age <- str_extract(race_result$horse_age, pattern = "\\d")
          race_result <- dplyr::mutate(race_result,horse_weight=as.character(str_extract_all(race_result$`horse_name/age`,"\\d{3}")))
          race_result <- dplyr::mutate(race_result,horse_weight_change=as.character(str_extract_all(race_result$`horse_name/age`,"\\([\\+|\\-]\\d+\\)|\\([\\d+]\\)")))
          race_result$horse_weight_change <- sapply(rm_round(race_result$horse_weight_change, extract=TRUE), paste, collapse="")
          race_result <- dplyr::mutate(race_result,brinker=as.character(str_extract_all(race_result$`horse_name/age`,"B")))
          race_result$brinker[race_result$brinker!="B"] <- "N"
          race_result <- race_result[-4]
          
          # ジョッキー
          race_result <- dplyr::mutate(race_result,jockey=as.character(str_extract_all(race_result$`jockey/weight`,"[ぁ-ん一-龠]+\\s[ぁ-ん一-龠]+|[:upper:].[ァ-ヶー]+")))
          race_result <- dplyr::mutate(race_result,jockey_weight=as.character(str_extract_all(race_result$`jockey/weight`,"\\d{2}")))
          race_result$jockey_weight_change <- 0
          race_result$jockey_weight_change[str_detect(race_result$`jockey/weight`,"☆")==1] <- 1
          race_result$jockey_weight_change[str_detect(race_result$`jockey/weight`,"△")==1] <- 2
          race_result$jockey_weight_change[str_detect(race_result$`jockey/weight`,"△")==1] <- 3
          race_result <- race_result[-4]
          
          # オッズと人気
          race_result <- dplyr::mutate(race_result,odds=as.character(str_extract_all(race_result$`popularity/odds`,"\\(.+\\)")))
          race_result <- dplyr::mutate(race_result,popularity=as.character(str_extract_all(race_result$`popularity/odds`,"\\d+[^(\\d+.\\d)]")))
          race_result$odds <- sapply(rm_round(race_result$odds, extract=TRUE), paste, collapse="")
          race_result <- race_result[-4]

次に、今取得したtable以外の情報も取り込むことにします。具体的には、レース名や天候、馬場状態、日付、競馬場などです。これらの情報はレース結果ページの上部に掲載されています。 f:id:osashimix:20180814215142p:plain

          # レース情報
          race_date <- race1 %>% html_nodes(xpath = "//div[@id = 'raceTitName']/p[@id = 'raceTitDay']") %>%
            html_text()
          race_name <- race1 %>% html_nodes(xpath = "//div[@id = 'raceTitName']/h1[@class = 'fntB']") %>%
            html_text()
          race_distance <- race1 %>% html_nodes(xpath = "//p[@id = 'raceTitMeta']") %>%
            html_text()
        
          race_result <- dplyr::mutate(race_result,race_date=as.character(str_extract_all(race_date,"\\d+年\\d+月\\d+日")))
          race_result$race_date <- str_replace_all(race_result$race_date,"年","/")
          race_result$race_date <- str_replace_all(race_result$race_date,"月","/")
          race_result$race_date <- as.Date(race_result$race_date)
          race_course <- as.character(str_extract_all(race_date,pattern = "札幌|函館|福島|新潟|東京|中山|中京|京都|阪神|小倉"))
          race_result$race_course <- race_course
          race_result <- dplyr::mutate(race_result,race_name=as.character(str_replace_all(race_name,"\\s","")))
          race_result <- dplyr::mutate(race_result,race_distance=as.character(str_extract_all(race_distance,"\\d+m")))
          race_type=as.character(str_extract_all(race_distance,pattern = "芝|ダート"))
          race_result$type <- race_type
          race_turn <- as.character(str_extract_all(race_distance,pattern = "右|左"))
          race_result$race_turn <- race_turn
          
          if(length(race1 %>% html_nodes(xpath = "//img[@class = 'spBg ryou']")) == 1){
            race_result$race_condition <- "良"
          } else if (length(race1 %>% 
                            html_nodes(xpath = "//img[@class = 'spBg yayaomo']")) == 1){
            race_result$race_condition <- "稍重"
          } else if (length(race1 %>%
                            html_nodes(xpath = "//img[@class = 'spBg omo']")) == 1){
            race_result$race_condition <- "重"
          } else if (length(race1 %>% 
                            html_nodes(xpath = "//img[@class = 'spBg furyou']")) == 1){
            race_result$race_condition <- "不良"
          } else race_result$race_condition <- "NA"
          
          if (length(race1 %>% html_nodes(xpath = "//img[@class = 'spBg hare']")) == 1){
            race_result$race_weather <- "晴れ"
          } else if (length(race1 %>% html_nodes(xpath = "//img[@class = 'spBg ame']")) == 1){
            race_result$race_weather <- "曇り"
          } else if (length(race1 %>% html_nodes(xpath = "//img[@class = 'spBg kumori']")) == 1){
            race_result$race_weather <- "雨"
          } else race_result$race_weather <- "その他"

次は各馬の情報です。 実はさきほど取得したtableの馬名はリンクになっており、そのリンクをたどると各馬の情報が取得できます（毛色や生年月日など）。

          horse_url <- race1 %>% html_nodes("a") %>% html_attr("href") 
          horse_url <- horse_url[str_detect(horse_url, pattern="directory/horse")==1] # 馬情報のリンクだけ抽出する
          
          for (l in 1:length(horse_url)){
            tryCatch(
              {
                horse1 <-  read_html(str_c("https://keiba.yahoo.co.jp",horse_url[l]))
                Sys.sleep(0.5)
                horse_name <- horse1 %>% html_nodes(xpath = "//div[@id = 'dirTitName']/h1[@class = 'fntB']") %>% 
                  html_text()
                horse <- horse1 %>% html_nodes(xpath = "//div[@id = 'dirTitName']/ul") %>% 
                  html_text()
                race_result$colour[race_result$horse_name==horse_name] <- as.character(str_extract_all(horse,"毛色：.+")) 
                race_result$owner[race_result$horse_name==horse_name] <- as.character(str_extract_all(horse,"馬主：.+"))
                race_result$farm[race_result$horse_name==horse_name] <- as.character(str_extract_all(horse,"生産者：.+"))
                race_result$locality[race_result$horse_name==horse_name] <- as.character(str_extract_all(horse,"産地：.+"))
                race_result$horse_birthday[race_result$horse_name==horse_name] <- as.character(str_extract_all(horse,"\\d+年\\d+月\\d+日"))
                race_result$father[race_result$horse_name==horse_name] <- horse1 %>% html_nodes(xpath = "//td[@class = 'bloodM'][@rowspan = '4']") %>% html_text()
                race_result$mother[race_result$horse_name==horse_name] <- horse1 %>% html_nodes(xpath = "//td[@class = 'bloodF'][@rowspan = '4']") %>% html_text()
              }
              , error = function(e){
                race_result$colour[race_result$horse_name==horse_name] <- NA 
                race_result$owner[race_result$horse_name==horse_name] <- NA
                race_result$farm[race_result$horse_name==horse_name] <- NA
                race_result$locality[race_result$horse_name==horse_name] <- NA
                race_result$horse_birthday[race_result$horse_name==horse_name] <- NA
                race_result$father[race_result$horse_name==horse_name] <- NA
                race_result$mother[race_result$horse_name==horse_name] <- NA
                }
            )
          }
          
          race_result$colour <- str_replace_all(race_result$colour,"毛色：","")
          race_result$owner <- str_replace_all(race_result$owner,"馬主：","")
          race_result$farm <- str_replace_all(race_result$farm,"生産者：","")
          race_result$locality <- str_replace_all(race_result$locality,"産地：","")
          #race_result$horse_birthday <- str_replace_all(race_result$horse_birthday,"年","/")
          #race_result$horse_birthday <- str_replace_all(race_result$horse_birthday,"月","/")
          #race_result$horse_birthday <- as.Date(race_result$horse_birthday)
          
          race_result <- dplyr::arrange(race_result,horse_number) # 馬番順に並べる

次にそのレースまでに獲得した賞金額を落としに行きます。これはレース結果のページの出馬表と書かれたリンクをたどるとアクセスできます。ここに賞金があるのでそれを取得します。

          yosou_url <- race1 %>% html_nodes("a") %>% html_attr("href") 
          yosou_url <- yosou_url[str_detect(yosou_url, pattern="denma")==1]
          
          if (length(yosou_url)==1){
          yosou1 <-  read_html(str_c("https://keiba.yahoo.co.jp",yosou_url)) 
          Sys.sleep(2)
          yosou <- yosou1 %>% html_nodes(xpath = "//td[@class = 'txC']") %>% as.character()
          prize <- yosou[grepl("万",yosou)==TRUE] %>% str_extract_all("\\d+万")
          prize <- t(do.call("data.frame",prize)) %>% as.character()
          race_result$prize <- prize
          race_result$prize <- str_replace_all(race_result$prize,"万","") %>% as.numeric()
          } else race_result$prize <- NA

取得した各レース結果を格納するdatasetというデータフレームを作成し、データを格納していきます。1年ごとにそれをSQLiteへ保存していきます。

          ## ファイル貯めるのかく
          if (k == 1 && i == 1 && j == 1){
            dataset <- race_result
          } else {
            dataset <- rbind(dataset,race_result)
          } # if文2の終わり
        }else
        {
          print("保存できませんでした") 
        }# if文1の終わり
      } # iループの終わり
    } # jループ終わり
  } # kループの終わり
  beep(3)
  write.csv(dataset,"race_result2.csv", row.names = FALSE)
  
  if (year == 1994){
    dbWriteTable(con, "race_result", dataset)
  } else {
    dbWriteTable(con, "temp", dataset)
    dbSendQuery(con, "INSERT INTO race_result select * from temp")
    dbSendQuery(con, "DROP TABLE temp")
  } # ifの終わり
} # yearループの終わり
end.time <- Sys.time()
print(str_c("処理時間は",end.time-start.time,"です。"))
beep(5)

options(warn = 1)

dbDisconnect(con)

以上です。結果はまた後程。

【日次GDP】Gianonne et. al. (2008)のマルチファクターモデルで四半期GDPを予想してみた

おはこんばんにちは。 前回、統計ダッシュボードからAPI接続で統計データを落とすという記事を投稿しました。

osashimix.hatenablog.com

今回はそのデータを、Gianonne et. al. (2008)のマルチファクターモデルにかけ、四半期GDPの予測を行いたいと思います。

Gianonne et. al. (2008)版マルチファクターモデル

元論文 http://dept.ku.edu/~empirics/Courses/Econ844/papers/Nowcasting%20GDP.pdf

前回の投稿でも書きましたが、この論文はGiannoneらが2008年にパブリッシュした論文です(JME)。彼らはアメリカの経済指標を用いて四半期GDPを日次で推計し、予測指標としての有用性を示しました。指標間の連動性(colinearity)を利用して、多数ある経済指標を2つのファクターに圧縮し、そのファクターを四半期GDPにフィッティングさせることによって高い予測性を実現しています。 まず、このモデルについてご紹介します。このモデルでは2段階推計を行います。まず主成分分析により経済統計を統計間の相関が0となるファクターへ変換します（

データセットの見える化・可視化といったらまずはこれ！～主成分分析(Principal Component Analysis, PCA)～

）。そして、その後の状態空間モデルでの推計で必要になるパラメータをOLS推計し、そのパラメータを使用してカルマンフィルタ＆カルマンスムーザーを回し、ファクターを推計しています。では、具体的な説明に移ります。 統計データをx_{i,t|v_j}と定義します。ここで、i=1,...,nは経済統計を表し（つまりnが全統計数）、t=1,...,T_{iv_j}は統計iのサンプル期間の時点を表しています（つまり、T_{iv_j}は統計iのその時点での最新データ日付を表す）。また、v_jはある時点j（2005年など）で得られる情報集合（vintage）を表しています。統計データx_{i,t|v_j}は以下のようにファクターf_{r,t}の線形結合で表すことができます（ここでrはファクターの数を表す）。

x_{i,t|v_j} = \mu_i + \lambda_{i1}f_{1,t} + ... + \lambda_{ir}f_{r,t} + \xi_{i,t|v_j} \tag{1}

\mu_iは定数項、 \lambda_{ir}はファクターローディング、 \xi_{i,t|v_j}はホワイトノイズの誤差項を表しています。これを行列形式で書くと以下のようになります。

x_{t|v_j} = \mu + \Lambda F_t + \xi_{t|v_j} = \mu + \chi_t + \xi_{t|v_j} \tag{2}

ここで、x_{t|v_j} = (x_{1,t|v_j}, ..., x_{n,t|v_j} )^{\mathrm{T}}、\xi_{t|v_j}=(\xi_{1,t|v_j}, ..., \xi_{n,t|v_j})^{\mathrm{T}}、 F_t = (f_{1,t}, ..., f_{r,t})^{\mathrm{T}}であり、\Lambdaは各要素が \lambda_{ij}のn\times r行列のファクターローディングを表しています。また、\chi_t = \Lambda F_tです。よって、ファクター F_tを推定するためには、データ x_{i,t|v_j}を以下のように基準化したうえで、分散共分散行列を計算し、その固有値問題を解けばよいという事になります。

\displaystyle z_{it} = \frac{1}{\hat{\sigma}_i}(x_{it} - \hat{\mu}_{it}) \tag{3}

ここで、\displaystyle \hat{\mu}_{it} = 1/T \sum_{t=1}^T x_{it}であり、 \hat{\sigma}_i = \sqrt{1/T \sum_{t=1}^T (x_{it}-\hat{\mu_{it}})^2}です（ここで Tはサンプル期間）。分散共分散行列 Sを以下のように定義します。

\displaystyle S = \frac{1}{T} \sum_{t=1}^T z_t z_t^{\mathrm{T}} \tag{4}

次に、 Sのうち、固有値を大きい順にr個取り出し、それを要素にした r \times r対角行列を D、それに対応する固有ベクトルを n \times r行列にしたものを Vと定義します。ファクター \tilde{F}_tは以下のように推計できます。

\tilde{F}_t = V^{\mathrm{T}} z_t \tag{5}

ファクターローディング \Lambdaと誤差項の共分散行列 \Psi = \mathop{\mathbb{E}} [\xi_t\xi^{\mathrm{T}}_t]は \tilde{F}_tを z_tに回帰することで推計します。

\displaystyle \hat{\Lambda} = \sum_{t=1}^T z_t \tilde{F}^{\mathrm{T}}_t (\sum_{t=1}^T\tilde{F}_t\tilde{F}^{\mathrm{T}}_t)^{-1} = V \tag{6}

\hat{\Psi} = diag(S - VDV) \tag{7}

注意して頂きたいのは、ここで推計した \tilde{F}_tは、以下の状態空間モデルでの推計に必要なパラメータを計算するための一時的な推計値であるという事です（２段階推計の１段階目という事）。

x_{t|v_j} = \mu + \Lambda F_t + \xi_{t|v_j} = \mu + \chi_t + \xi_{t|v_j} \tag{2}

F_t = AF_{t-1} + u_t \tag{8}

ここで、u_tは平均0、分散Hのホワイトノイズです。再掲している(2)式が観測方程式、(8)式が遷移方程式となっています。推定すべきパラメータは \Lambda、 \Psi以外にAとHがあります（ \mu = 0としています）。Aは主成分分析により計算した \tilde{F}_tをVAR(1)にかけることで推定します。

\hat{A} = \sum_{t=2}^T\tilde{F}_t\tilde{F}^{\mathrm{T}}_{t-1} (\sum_{t=2}^T\tilde{F}_{t-1}\tilde{F}^{\mathrm{T}}_{t-1})^{-1} \tag{9}

Hは今推計したVAR(1)の誤差項の共分散行列から計算します。これで必要なパラメータの推定が終わりました。次にカルマンフィルタを回します。カルマンフィルタに関しては

シンプルなモデルとイラストでカルマンフィルタを直観的に理解してみる

を参考にしてください。わかりやすいです。これで最終的に\hat{F}_{t|v_j}の推計ができるわけです。 GDPがこれらのファクターで説明可能であり（つまり固有の変動がない）、GDPと月次経済指標がjointly normalであれば以下のような単純なOLS推計でGDPを予測することができます。もちろん月次経済指標の方が早く公表されるので、内生性の問題はないと考えられます。

\hat{y}_{3k|v_j} = \alpha + \beta^{\mathrm{T}} \hat{F}_{3k|v_j} \tag{10}

ここで、3kは四半期の最終月を示しています（3月、6月など）。 \hat{y}_{3k|v_j}はj時点で得られる情報集合 v_jでの四半期GDPを表しており、 \hat{F}_{3k|v_j}はその時点で推定したファクターを表しています（四半期最終月の値だけを使用している点に注意）。これで推計方法の説明は終わりです。

Rで実装する

では実装します。前回記事で得られたデータ（dataset1）が読み込まれている状態からスタートします。まず、主成分分析でファクターを計算します。なお、前回の記事で3ファクターの累積寄与度が80%を超えたため、今回もファクター数は3にしています。

#------------------------
# Giannone et. al. 2008 
#------------------------

# 主成分分析でファクターを計算
f <- 3 # ファクター数を定義
z <- scale(dataset1) # (1)式のデータセット基準化
for (i in 1:nrow(z)){
  eval(parse(text = paste("S_i <- z[i,]%*%t(z[i,])",sep = "")))
  if (i==1){
    S <- S_i
  }else{
    S <- S + S_i
  }
}
S <- (1/nrow(z))*S # 分散共分散行列を計算 (4)式
gamma <- eigen(S) 
D <- diag(gamma$values[1:f])
V <- gamma$vectors[,1:f]
F_t <- matrix(0,nrow(z),f)
for (i in 1:nrow(z)){
  eval(parse(text = paste("F_t[",i,",]<- z[",i,",]%*%V",sep = ""))) # (5)式を実行
}
F_t.xts <- xts(F_t,order.by = as.Date(row.names(z)))
plot.zoo(F_t.xts,col = c("red","blue","green","yellow","purple"),plot.type = "single") # 時系列プロット
lambda_hat <- V
psi <- diag(S-V%*%D%*%t(V)) # (7)式
R <- diag(diag(cov(z-z%*%V%*%t(V))))

推計したファクター \tilde{F}_tの時系列プロットは以下のようになり、前回princomp関数で計算したファクターと完全一致します（じゃあprincompでいいやんと思われるかもしれませんが実装しないと勉強になりませんので）。 f:id:osashimix:20180716144839p:plain 次に、VAR(1)を推計し、パラメータを取り出します。

# VAR(1)モデルを推計
a <- matrix(0,f,f)
b <- matrix(0,f,f)
for(t in 2:nrow(z)){
  a <- a + F_t[t,]%*%t(F_t[t-1,])
  b <- b + F_t[t-1,]%*%t(F_t[t-1,])
}
b_inv <- solve(b)
A_hat <- a%*%b_inv # (9)式
> A_hat
           [,1]        [,2]        [,3]
[1,] 0.99477825  0.01696902 -0.04746482
[2,] 0.00456925  0.99227689  0.03772908
[3,] 0.00596195 -0.02083938  0.97124953
e <- numeric(f)
for (t in 2:nrow(F_t)){
  e <- e + F_t[t,]-F_t[t-1,]%*%A_hat
}
H <- t(e)%*%e # 残差の分散共分散行列を計算
Q <- diag(1,f,f)
Q[1:f,1:f] <- H
> H
          [,1]     [,2]      [,3]
[1,] 24.489201 32.11384  9.335003
[2,] 32.113841 42.11239 12.241429
[3,]  9.335003 12.24143  3.558397

VAR(1)に関してもvar関数とパラメータの数値が一致することを確認済みです。いよいよカルマンフィルタを実行します。

# カルマンフィルタを実行
a <- which(dataset$publication == "2005-04-01") # サンプル開始期間を2005年に設定。
dataset2 <- dataset[a:nrow(dataset),]
rownames(dataset2) <- dataset2$publication
dataset2 <- dataset2[,-2]
z <- scale(dataset2) # zは基準化されたサンプルデータ
RR <- array(0,dim = c(ncol(z),ncol(z),nrow(z))) # RRは観測値の分散行列（相関はないと仮定）
for(i in 1:nrow(z)){
  miss <- is.na(z[i,])
  R_temp <- diag(R)
  R_temp[miss] <- 1e+32 # 欠損値の分散は無限大にする
  RR[,,i] <- diag(R_temp)
}
zz <- z; zz[is.na(z)] <- 0 # 欠損値（NA）に0を代入（計算結果にはほとんど影響しない）。
a_t <- matrix(0,nrow(zz),f) # a_tは状態変数の予測値
a_tt <- matrix(0,nrow(zz),f) # a_ttは状態変数の更新後の値
a_tt[1,] <- F_t[1,] # 状態変数の初期値には主成分分析で推計したファクターを使用
sigma_t <- array(0,dim = c(f,f,nrow(zz))) # sigma_tは状態変数の分散の予測値
sigma_tt <- array(0,dim = c(f,f,nrow(zz))) # sigma_tは状態変数の分散の更新値
p <- ginv(diag(nrow(kronecker(A_hat,A_hat)))-kronecker(A_hat,A_hat))
sigma_tt[,,1] <- matrix(p,3,3) # 状態変数の分散の初期値はVAR(1)の推計値から計算
y_t <- matrix(0,nrow(zz),ncol(zz)) # y_tは観測値の予測値
K_t <- array(0,dim = c(f,ncol(zz),nrow(zz))) # K_tはカルマンゲイン
data.m <- as.matrix(dataset2)
# カルマンフィルタを実行
for (t in 2:nrow(zz)){
  a_t[t,] <- A_hat%*%a_tt[t-1,]
  sigma_t[,,t] <- A_hat%*%sigma_tt[,,t-1]%*%t(A_hat) + Q
  y_t[t,] <- V%*%a_t[t,]
  S_t <- V%*%sigma_tt[,,t-1]%*%t(V)+RR[,,t]
  GG <- t(V)%*%diag(1/diag(RR[,,t]))%*%V
  Sinv <- diag(1/diag(RR[,,t])) - diag(1/diag(RR[,,t]))%*%V%*%ginv(diag(nrow(A_hat))+sigma_t[,,t]%*%GG)%*%sigma_t[,,t]%*%t(V)%*%diag(1/diag(RR[,,t]))
  K_t[,,t] <- sigma_t[,,t]%*%t(V)%*%Sinv
  a_tt[t,] <- a_t[t,] + K_t[,,t]%*%(zz[t,]-y_t[t,])
  sigma_tt[,,t] <- sigma_t[,,t] - K_t[,,t]%*%V%*%sigma_tt[,,t-1]%*%t(V)%*%t(K_t[,,t])
  }

F.xts <- xts(a_tt,order.by = as.Date(rownames(data.m)))
plot.zoo(F.xts, col = c("red","blue","green","yellow","purple"),plot.type = "single") # 得られた推計値を時系列プロット

カルマンフィルタにより推計したファクターの時系列プロットが以下です。遷移方程式がAR(1)だったからかかなり平準化された値となっています。 f:id:osashimix:20180716151958p:plain では、この得られたファクターをOLSにかけます。

# 得られたファクターとGDPをOLSにかける
F_q <- as.data.frame(a_tt[seq(3,nrow(a_tt),3),]) # 四半期の終わり月の値だけを引っ張ってくる 
colnames(F_q) <- c("factor1","factor2","factor3")
colnames(GDP) <- c("publication","GDP")
t <- which(GDP$publication=="2005-04-01")
t2 <- which(GDP$publication=="2015-01-01") # 2005-2q~2015-1qまでのデータが学習データ、それ以降がテストデータ
GDP_q <- GDP[t:nrow(GDP),]
dataset.q <- cbind(GDP_q[1:(t2-t),],F_q[1:(t2-t),])
test <- lm(GDP~factor1 + factor2 + factor3,data=dataset.q)
summary(test)
> summary(test)

Call:
lm(formula = GDP ~ factor1 + factor2 + factor3, data = dataset.q)

Residuals:
    Min      1Q  Median      3Q     Max 
-6032.0  -485.8  -139.5   920.4  3367.5 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 127388.6      772.0 165.019  < 2e-16 ***
factor1      -1530.2      245.2  -6.242 7.11e-07 ***
factor2      -1425.0      316.0  -4.509 9.28e-05 ***
factor3        238.0      287.6   0.828    0.414    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1716 on 30 degrees of freedom
Multiple R-squared:  0.6359, Adjusted R-squared:  0.5995 
F-statistic: 17.47 on 3 and 30 DF,  p-value: 9.51e-07

out_of_sample <- cbind(GDP_q[(t2-t+1):nrow(GDP_q),],F_q[(t2-t+1):nrow(GDP_q),]) # out of sampleのデータセットを作成
test.pred <-  predict(test, out_of_sample, interval="prediction")
pred.GDP.xts <- xts(cbind(test.pred[,1],out_of_sample$GDP),order.by = out_of_sample$publication)
plot.zoo(pred.GDP.xts,col = c("red","blue"),plot.type = "single") # 予測値と実績値を時系列プロット

OLSの推計結果はfactor1（赤）とfactor2（青）が有意との結果。前回の投稿でも言及したように、factor1（赤）はリスクセンチメントを表していそうなので、係数の符号が負であることは頷ける。ただし、factor2（青）も符号が負なのではなぜなのか…。このファクターは生産年齢人口など経済の潜在能力を表していると思っていたのに。かなり謎。まあとりあえず予測に移りましょう。このモデルを使用したGDPの予測値と実績値の推移はいかのようになりました。直近の精度は悪くない？ f:id:osashimix:20180716152708p:plain というか、これ完全に単位根の問題を無視してOLSしてしまっているな。ファクターもGDPも完全に単位根を持つけど念のため単位根検定をかけてみます。

> adf.test(F_q$factor1)

    Augmented Dickey-Fuller Test

data:  F_q$factor1
Dickey-Fuller = -1.4532, Lag order = 3, p-value = 0.7946
alternative hypothesis: stationary
> adf.test(F_q$factor2)

    Augmented Dickey-Fuller Test

data:  F_q$factor2
Dickey-Fuller = -2.388, Lag order = 3, p-value = 0.4183
alternative hypothesis: stationary

> adf.test(F_q$factor3)

    Augmented Dickey-Fuller Test

data:  F_q$factor3
Dickey-Fuller = -2.728, Lag order = 3, p-value = 0.2815
alternative hypothesis: stationary

> adf.test(GDP_q$GDP)

    Augmented Dickey-Fuller Test

data:  GDP_q$GDP
Dickey-Fuller = -2.0259, Lag order = 3, p-value = 0.564
alternative hypothesis: stationary

はい。全部単位根もってました…。階差をとったのち、単位根検定を行います。

GDP <- GDP %>% mutate(growth.rate=(GDP/lag(GDP)-1)*100)
F_q <- F_q %>% mutate(f1.growth.rate=(factor1/lag(factor1)-1)*100,
                      f2.growth.rate=(factor2/lag(factor2)-1)*100,
                      f3.growth.rate=(factor3/lag(factor3)-1)*100)
> adf.test(GDP$growth.rate[2:NROW(GDP$growth.rate)])

    Augmented Dickey-Fuller Test

data:  GDP$growth.rate[2:NROW(GDP$growth.rate)]
Dickey-Fuller = -4.6703, Lag order = 5, p-value = 0.01
alternative hypothesis: stationary
> adf.test(F_q$f1.growth.rate[2:NROW(F_q$f1.growth.rate)])

    Augmented Dickey-Fuller Test

data:  F_q$f1.growth.rate[2:NROW(F_q$f1.growth.rate)]
Dickey-Fuller = -3.4292, Lag order = 3, p-value = 0.06092
alternative hypothesis: stationary

> adf.test(F_q$f2.growth.rate[2:NROW(F_q$f2.growth.rate)])

    Augmented Dickey-Fuller Test

data:  F_q$f2.growth.rate[2:NROW(F_q$f2.growth.rate)]
Dickey-Fuller = -4.0271, Lag order = 3, p-value = 0.01518
alternative hypothesis: stationary

> adf.test(F_q$f3.growth.rate[2:NROW(F_q$f3.growth.rate)])

    Augmented Dickey-Fuller Test

data:  F_q$f3.growth.rate[2:NROW(F_q$f3.growth.rate)]
Dickey-Fuller = -3.7663, Lag order = 3, p-value = 0.02769
alternative hypothesis: stationary

factor1だけは5%有意水準で帰無仮説を棄却できない…。困りました。有意水準を10%ということにして、とりあえず階差でOLSしてみます。

test1 <- lm(growth.rate~f1.growth.rate + f2.growth.rate + f3.growth.rate,data=dataset.q)
summary(test1)
>summary(test1)
Call:
lm(formula = growth.rate ~ f1.growth.rate + f2.growth.rate + 
    f3.growth.rate, data = dataset.q)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6702 -0.5511  0.0774  0.6978  2.5696 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)    -0.015026   0.195777  -0.077   0.9393  
f1.growth.rate -0.001975   0.001344  -1.470   0.1508  
f2.growth.rate  0.007496   0.004568   1.641   0.1100  
f3.growth.rate  0.003282   0.001371   2.394   0.0223 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.163 on 34 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.233,  Adjusted R-squared:  0.1653 
F-statistic: 3.443 on 3 and 34 DF,  p-value: 0.02739

推計結果がわるくなりました…。予測値を計算し、実績値とプロットしてみます。 f:id:osashimix:20180716162709p:plain ん～、これはやり直しですね。今日はここまでで勘弁してください…。

【日次GDP】日次GDP推計に使用する経済統計を統計ダッシュボードから集めてみた。

日次GDP

おはこんばんちわ。

最近、競馬ばっかりやってましたが、そろそろ本業のマクロの方もやらないとなということで今回は日次GDP推計に使用するデータを総務省が公開している統計ダッシュボードから取ってきました。 そもそも、前の記事では四半期GDP速報の精度が低いことをモチベーションに高頻度データを用いてより精度の高い予測値をはじき出すモデルを作れないかというテーマで研究を進めていました。しかし、先行研究を進めていくうちに、どうやら大規模な経済指標を利用することで日次で四半期GDPの予測値を計算することが可能であることが判明しました。しかも、精度も良い(米国ですが)ということで、なんとかこの方向で研究を進めていけないかということになりました。

先行研究と具体的にやりたいこと

先行研究

http://dept.ku.edu/~empirics/Courses/Econ844/papers/Nowcasting%20GDP.pdf

Giannoneらが2008年にパブリッシュした論文です(JME)。彼らはアメリカの経済指標を用いて四半期GDPを日次で推計し、予測指標としての有用性を示しました。指標間の連動性(colinearity)を利用して、多数ある経済指標をいくつかのファクターに圧縮し、そのファクターを四半期GDPにフィッティングさせることによって高い予測性を実現しました。なお、ファクターの計算にはカルマンスムージングを用いています(詳しい推計方法は論文&Technical Appendixを参照)。理論的な定式化は無いのですが、なかなか当たります。そもそも私がこの研究に興味を持ったのは、以下の本を立ち読みした際に参考文献として出てきたからで、いよいよ運用機関などでも使用され始めるのかと思い、やっておこうと思った次第です。

実践 金融データサイエンス 隠れた構造をあぶり出す6つのアプローチ

実践 金融データサイエンス 隠れた構造をあぶり出す6つのアプローチ

とりあえずはGiannoneの日本版をやろうかなと思っています。実はこの後に、ファクターモデルとDSGEを組み合わせたモデルがありましてそこまで発展させたいなーなんて思っておりますが。とにかく、ファクターを計算するための経済統計が必要ですので、今回はそれを集めてきたというのがこの記事の趣旨です。

統計ダッシュボードからのデータの収集

政府や日銀が公表しているデータの一部は統計ダッシュボードから落とすことができます。これは総務省統計局が提供しているもので、これまで利用しにくかった経済統計をより身近に使用してもらおうというのが一応のコンセプトとなっています。似たものに総務省統計局が提供しているestatがありますが、日銀の公表データがなかったり、メールアドレスの登録が必要だったりと非常に使い勝手が悪いです(個人的感想)。ただ、estatにはestatapiというRパッケージがあり、データを整形するのは比較的容易であると言えます。今回、統計ダッシュボードを選択した理由はそうは言っても日銀のデータがないのはダメだろうという理由で、データの整形に関しては関数を組みました。 そもそも統計ダッシュボードは経済統計をグラフなどで見て楽しむ？ものですが、私のような研究をしたい者を対象にAPIを提供してくれています。取得できるデータは大きく分けて6つあります。 f:id:osashimix:20180714160029p:plain

やり方は簡単で、ベースのurlと欲しい統計のIDをGET関数で渡すことによって、データを取得することができます。公式にも以下のように書かれています。

基本的な使い方としては、まず①「統計メタ情報（系列）取得」で取得したいデータの[系列コード]を検索し、 その後⑥「統計データ取得」で[系列コード]を検索条件に指定し、その系列の情報を取得します。 （②③④⑤は補助的な情報として独立して取得できるようにしています。データのみ必要な場合は当該機能は不要です。） 具体的な使い方は、以下の「WebAPIの詳細仕様」に記載する[ベースURL]に検索条件となる[パラメータ]を"&"で連結し、HTTPリクエスト（GET）を送信することで目的のデータを取得できます。 各パラメータは「パラメータ名=値」のように名称と値を'='で結合し、複数のパラメータを指定する場合は「パラメータ名=値&パラメータ名=値&…」のようにそれぞれのパラメータ指定を'&'で結合してください。 また、パラメータ値は必ずURLエンコード(文字コードUTF-8)してから結合してください。

今回も以下の文献を参考にデータを取ってきたいと思います。

Rによるスクレイピング入門

Rによるスクレイピング入門

まず、最初にこのAPIからデータを取得し、得られた結果を分析しやすいように整形する関数を定義したいと思います。

library(httr)
library(dplyr)
library(stringr)
library(xts)
library(ggplot2)
library(seasonal)

# 関数を定義
get_dashboard <- function(ID){
  base_url <- "https://dashboard.e-stat.go.jp/api/1.0/JsonStat/getData?"
  res <- GET(
    url = base_url,
    query = list(
      IndicatorCode=ID
    )
  )
  result <- content(res)
  x <- result$link$item[[1]]$value
  x <- t(do.call("data.frame",x))
  date_x <- result$link$item[[1]]$dimension$Time$category$label
  date_x <- t(do.call("data.frame",date_x))
  date_x <- str_replace_all(date_x, pattern="年", replacement="/")
  date_x <- str_replace_all(date_x, pattern="月", replacement="")
  date_x <- as.Date(gsub("([0-9]+)/([0-9]+)", "\\1/\\2/1", date_x))
  date_x <- as.Date(date_x, format = "%m/%d/%Y")
  date_x <- as.numeric(date_x)
  date_x <- as.Date(date_x, origin="1970-01-01")
  #x <- cbind(x,date_x)
  x <- data.frame(x)
  x[,1] <- as.character(x[,1])%>%as.numeric(x[,1])
  colnames(x) <- c(result$link$item[[1]]$label)
  x <- x %>% mutate("publication" = date_x)
  return(x)
}

まずベースのurlを定義しています。今回はデータが欲しいので⑥統計データのベースurlを使用します（ 統計ダッシュボード - API）。次にベースurlと統計ID（IndicatorCode）をGET関数で渡し、結果を取得しています。統計IDについてはエクセルファイルで公開されています。得られた結果の中身（リスト形式）をresultに格納し、リストの深層にある原数値データ（value）をxに格納します。原数値データもリスト形式なので、それをdo.callでデータフレームに変換しています。次に、データ日付を取得します。resultの中を深くたどるとTime$category$labelというデータがあり、そこに日付データが保存されているので、それをdate_xに格納し、同じようにデータフレームへ変換します。データの仕様上、日付は「yyyy年mm月」になっていますが、これだとRは日付データとして読み取ってくれないので、str_replace_all等で変換したのち、Date型に変換しています。列名にデータ名（result$link$item1$label）をつけ、データ日付をxに追加したら完成です。 そのほか、data_connectという関数も定義しています。これはデータ系列によれば、たとえば推計方法の変更などで1980年～2005年の系列と2003年～2018年までの系列の2系列があるようなデータも存在し、この2系列を接続するための関数です。これは単純に接続しているだけなので、説明は省略します。

data_connect <- function(x){
  a <- min(which(x[,ncol(x)] != "NA"))
  b <- x[a,ncol(x)]/x[a,1]
  c <- x[1:a-1,1]*b
return(c)
}

では、実際にデータを取得していきます。今回取得するデータは月次データとなっています。これは統計dashboardが月次以下のデータがとれないからです。なので、例えば日経平均などは月末の終値を引っ張っています。ただし、GDPは四半期データとなっています。さきほど定義したget_dashboardの使用方法は簡単で、引数に統計ダッシュボードで公開されている統計IDを入力するだけでデータが取れます。今回使用するデータを以下の表にまとめました。 f:id:osashimix:20180714212330p:plain

# データを取得
Nikkei <- get_dashboard("0702020501000010010")
callrate <- get_dashboard("0702020300000010010")
kikai <- get_dashboard("0701030000000010010")
kigyo.bukka <- get_dashboard("0703040300000090010")
money.stock1 <- get_dashboard("0702010201000010030")
money.stock2 <- get_dashboard("0702010202000010030")
money.stock <- dplyr::full_join(money.stock1,money.stock2,by="publication")
c <- data_connect(money.stock)
a <- min(which(money.stock[,ncol(money.stock)] != "NA"))
money.stock[1:a-1,ncol(money.stock)] <- c
money.stock <- money.stock[,c(2,3)]
cpi <- get_dashboard("0703010401010090010")
export.price <- get_dashboard("0703050301000090010")
import.price <- get_dashboard("0703060301000090010")
import.price$`輸出物価指数（総平均）（円ベース）2015年基準` <- NULL
public.expenditure1 <- get_dashboard("0802020200000010010")
public.expenditure2 <- get_dashboard("0802020201000010010")
public.expenditure <- dplyr::full_join(public.expenditure1,public.expenditure2,by="publication")
c <- data_connect(public.expenditure)
a <- min(which(public.expenditure[,ncol(public.expenditure)] != "NA"))
public.expenditure[1:a-1,ncol(public.expenditure)] <- c
public.expenditure <- public.expenditure[,c(2,3)]
export.service <- get_dashboard("1601010101000010010")
working.population <- get_dashboard("0201010010000010020")
unemployment_rate <- get_dashboard("0301010000020020010")
yukoukyuujinn <- get_dashboard("0301020001000010010")
hours_worked <- get_dashboard("0302010000000010000")
nominal.wage <- get_dashboard("0302020000000010000") 
iip <- get_dashboard("0502070101000090010")
shukka.shisu <- get_dashboard("0502070102000090010")
zaiko.shisu <- get_dashboard("0502070103000090010")
sanji.sangyo <- get_dashboard("0603100100000090010")
retail.sells <- get_dashboard("0601010201010010000")
GDP1 <- get_dashboard("0705020101000010000")
GDP2 <- get_dashboard("0705020301000010000")
GDP <- dplyr::full_join(GDP1,GDP2,by="publication")
c <- data_connect(GDP)
a <- min(which(GDP[,ncol(GDP)] != "NA"))
GDP[1:a-1,ncol(GDP)] <- c
GDP <- GDP[,c(2,3)]
yen <- get_dashboard("0702020401000010010")

今取得したデータは原数値系列のデータが多いので、それらは季節調整をかけます。なぜ季節調整済みのデータを取得しないのかというとそれらのデータは何故か極端にサンプル期間が短くなってしまうからです。ここらへんは使い勝手が悪いです。

# 季節調整をかける
Sys.setenv(X13_PATH = "C:/winx13_V2.5/WinX13/x13as")
checkX13()
# 季節調整済み系列を抜き出す関数を作成
seasoning <- function(data,i,start.y,start.m){
  timeseries <- ts(data[,i],frequency = 12,start=c(start.y,start.m))
  m <- seas(timeseries)
  return(m$series$s11)
}
k <- seasoning(kikai,1,2005,4)
kikai$`機械受注額（船舶・電力を除く民需）` <- as.numeric(k)
k <- seasoning(kigyo.bukka,1,1960,1)
kigyo.bukka$`国内企業物価指数（総平均）2015年基準` <- as.numeric(k)
k <- seasoning(cpi,1,1970,1)
cpi$`消費者物価指数（生鮮食品を除く総合）2015基準` <- as.numeric(k)
k <- seasoning(export.price,1,1960,1)
export.price$`輸出物価指数（総平均）（円ベース）2015年基準` <- as.numeric(k)
k <- seasoning(import.price,1,1960,1)
import.price$`輸入物価指数（総平均）（円ベース）2015年基準` <- as.numeric(k)
k <- seasoning(public.expenditure,2,2004,4)
public.expenditure$公共工事受注額 <- as.numeric(k)
k <- seasoning(export.service,1,1996,1)
export.service$`貿易・サービス収支` <- as.numeric(k)
k <- seasoning(unemployment_rate,1,1968,1)
unemployment_rate$`完全失業率（男女計）` <- as.numeric(k)
k <- seasoning(yukoukyuujinn,1,1963,1)
yukoukyuujinn$有効求人倍率 <- as.numeric(k)
k <- seasoning(hours_worked,1,1990,1)
hours_worked$総実労働時間 <- as.numeric(k)
k <- seasoning(nominal.wage,1,1990,1)
nominal.wage$現金給与総額 <- as.numeric(k)
k <- seasoning(iip,1,1978,1)
iip$`鉱工業生産指数　2010年基準` <- as.numeric(k)
k <- seasoning(shukka.shisu,1,1990,1)
shukka.shisu$`鉱工業出荷指数　2010年基準` <- as.numeric(k)
k <- seasoning(zaiko.shisu,1,1990,1)
zaiko.shisu$`鉱工業在庫指数　2010年基準` <- as.numeric(k)
k <- seasoning(sanji.sangyo,1,1988,1)
sanji.sangyo$`第３次産業活動指数　2010年基準` <- as.numeric(k)
k <- seasoning(retail.sells,1,1980,1)
retail.sells$`小売業販売額（名目）` <- as.numeric(k)
k <- seasoning(household.consumption,1,2010,1)
household.consumption$`二人以上の世帯　消費支出（除く住居等）` <- as.numeric(k)
GDP.ts <- ts(GDP[,2],frequency = 4,start=c(1980,1))
m <- seas(GDP.ts)
GDP$`国内総生産（支出側）（実質）2011年基準` <- m$series$s11

ここでは詳しく季節調整のかけ方は説明しません。x13arimaを使用しています。上述のコードを回す際はx13arimaがインストールされている必要があります。以下の記事を参考にしてください。

sinhrks.hatenablog.com

では、データ日付を基準に落としてきたデータを結合し、データセットを作成します。

# データセットに結合
dataset <- dplyr::full_join(kigyo.bukka,callrate,by="publication")
dataset <- dplyr::full_join(dataset,kikai,by="publication")
dataset <- dplyr::full_join(dataset,Nikkei,by="publication")
dataset <- dplyr::full_join(dataset,money.stock,by="publication")
dataset <- dplyr::full_join(dataset,cpi,by="publication")
dataset <- dplyr::full_join(dataset,export.price,by="publication")
dataset <- dplyr::full_join(dataset,import.price,by="publication")
dataset <- dplyr::full_join(dataset,public.expenditure,by="publication")
dataset <- dplyr::full_join(dataset,export.service,by="publication")
dataset <- dplyr::full_join(dataset,working.population,by="publication")
dataset <- dplyr::full_join(dataset,unemployment_rate,by="publication")
dataset <- dplyr::full_join(dataset,yukoukyuujinn,by="publication")
dataset <- dplyr::full_join(dataset,hours_worked,by="publication")
dataset <- dplyr::full_join(dataset,nominal.wage,by="publication")
dataset <- dplyr::full_join(dataset,iip,by="publication")
dataset <- dplyr::full_join(dataset,shukka.shisu,by="publication")
dataset <- dplyr::full_join(dataset,zaiko.shisu,by="publication")
dataset <- dplyr::full_join(dataset,sanji.sangyo,by="publication")
dataset <- dplyr::full_join(dataset,retail.sells,by="publication")
dataset <- dplyr::full_join(dataset,yen,by="publication")
colnames(dataset) <- c("DCGPI","publication","callrate","Machinery_Orders",
                       "Nikkei225","money_stock","CPI","export_price",
                       "import_price","public_works_order",
                       "trade_service","working_population",
                       "unemployment_rate","active_opening_ratio","hours_worked",
                       "wage","iip_production","iip_shipment","iip_inventory",
                       "ITIA","retail_sales","yen")

最後に列名をつけています。datasetはそれぞれのデータの公表開始時期が異なるために大量のNAを含むデータフレームとなっているので、NAを削除するために最もデータの開始時期が遅い機械受注統計に合わせてデータセットを再構築します。

a <- min(which(dataset$Machinery_Orders != "NA"))
dataset1 <- dataset[a:nrow(dataset),]
dataset1 <- na.omit(dataset1)
rownames(dataset1) <- dataset1$publication
dataset1 <- dataset1[,-2]
dataset1.xts <- xts(dataset1,order.by = as.Date(rownames(dataset1)))

各データの推移は以下のようになっています（見にくくてすみません）。 f:id:osashimix:20180714221953p:plain f:id:osashimix:20180714223037p:plain f:id:osashimix:20180714223107p:plain f:id:osashimix:20180714222734p:plain これでとりあえずデータの収集は終わりました。

得られたデータを主成分分析にかけてみる

本格的な分析はまた今後にしたいのですが、データを集めるだけでは面白くないので、Gianonneらのように主成分分析を行いたいと思います。主成分分析をこれまでに学んだことのない方は以下を参考にしてください。個人的にはわかりやすいと思っています。

datachemeng.com

では主成分分析を実行してみます。Rではprincomp関数を使用することで非常に簡単に主成分分析を行うことができます。

# 主成分分析を実行
factor.pca <- princomp(~.,cor = TRUE,data = dataset1) # cor = TRUEでデータの基準化を自動で行ってくれる。
summary(factor.pca)
> summary(factor.pca)
Importance of components:
                          Comp.1    Comp.2    Comp.3    Comp.4     Comp.5     Comp.6
Standard deviation     2.7811619 2.6222197 1.5406503 1.0303033 0.94327309 0.83277036
Proportion of Variance 0.3683267 0.3274303 0.1130287 0.0505488 0.04236972 0.03302412
Cumulative Proportion  0.3683267 0.6957570 0.8087857 0.8593346 0.90170427 0.93472839
                           Comp.7     Comp.8    Comp.9     Comp.10     Comp.11
Standard deviation     0.64851803 0.55413673 0.4584591 0.336977547 0.326749642
Proportion of Variance 0.02002741 0.01462226 0.0100088 0.005407327 0.005084063
Cumulative Proportion  0.95475580 0.96937806 0.9793869 0.984794189 0.989878252
                           Comp.12    Comp.13    Comp.14     Comp.15      Comp.16
Standard deviation     0.249285611 0.20512903 0.18200901 0.152523504 0.1365294632
Proportion of Variance 0.002959206 0.00200371 0.00157749 0.001107782 0.0008876331
Cumulative Proportion  0.992837458 0.99484117 0.99641866 0.997526440 0.9984140726
                            Comp.17      Comp.18      Comp.19      Comp.20
Standard deviation     0.1167617026 0.0824373195 0.0787936993 0.0651688158
Proportion of Variance 0.0006492045 0.0003236148 0.0002956403 0.0002022369
Cumulative Proportion  0.9990632772 0.9993868920 0.9996825323 0.9998847692
                            Comp.21
Standard deviation     0.0491919338
Proportion of Variance 0.0001152308
Cumulative Proportion  1.0000000000

screeplot(factor.pca)
pc <- predict(factor.pca,dataset1)[,1:3] # 主成分を計算
pc.xts <- xts(pc,order.by = as.Date(rownames(dataset1)))
plot.zoo(pc.xts,col=c("red","blue","green","purple","yellow"),plot.type = "single") # 主成分を時系列プロット

第3主成分まででデータの約80％が説明できる結果を得たので、第3主成分までのプロットをお見せします。第1主成分（赤）はリーマンショックや東日本大震災、消費税増税のあたりで急上昇しています。ゆえに経済全体のリスクセンチメントを表しているのではないかと思っています。第2主成分（青）と第3主成分（緑）はリーマンショックのあたりで大きく落ち込んでいることは共通していますが2015年～現在の動きが大きく異なっています。また、第2主成分（青）はサンプル期間を通して過去トレンドを持つことから日本経済の潜在能力のようなものを表しているのではないでしょうか（そうするとリーマンショックまで上昇傾向にあることが疑問なのですが）。第3主成分（緑）はいまだ解読不能です（物価＆為替動向を表しているのではないかと思っています）。とりあえず今日はこれまで。次回はGianonne et. al.(2008)の日本版の再現を行いたいと思います。

f:id:osashimix:20180714224518p:plain