雇用統計予測モデルの構築


## 概要
ご存知の通り、米国の雇用統計はFXの世界では最重要指標として注目されます。基本的に毎月第一金曜日の日本時間22:00くらいに発表され、その結果によって為替レートが大きく動きます。
一般に雇用統計では民間予測と発表値の乖離がどのくらいあるかが重要になってきます。ですので、より正しい予測をすることは重要です。

今回は色々な指標から発表前に予測が可能な雇用統計予測モデルを構築したいと思います。
モデル構築の注意点として、使用する指標が発表前に入手可能か、またどのくらい前に入手可能かを気にする必要があります。
十分注意しているつもりですが、もし無理じゃね？という指標があればご指摘願えると幸いです。

## 参考
なお、以下の記事にインスパイアされてます。

http://fx-kirin.com/python/nfp-linear-reggression-model/

## データ

上記の参考記事では結構たくさんの指標が使われています。

### 説明変数候補

* ニューヨーク連銀製造業景気指数の「雇用」
* フィラデルフィア連銀製造業景況指数の「雇用」
* ミシガン大学消費者信頼感指数
* カンザスシティ連銀製造業景況指数の「雇用」
* ダラス連銀製造業活動指数の「雇用」
* リッチモンド連銀製造業景況指数の「雇用」
* 消費者信頼感指数の「雇用不十分」
* 消費者信頼感指数の「雇用困難」
* 米ドル/円変化率（１カ月前）
* ＩＳＭ製造業景況指数の「雇用」
* Ｓ＆Ｐ500指数変化率（１カ月前）
* ＡＤＰ雇用統計
* チャレンジャー人員削減予定数

このうち

* ニューヨーク連銀製造業景気指数の「雇用」
* リッチモンド連銀製造業景況指数の「雇用」
* ＡＤＰ雇用統計
* Dow Jones月次変化率 (SPの代わり)

あたりを使用します。
加えて独自指標を何個か使います。

### 被説明変数

被説明変数は雇用統計の非農業部門新規雇用者数前月差を使用します。

URL : http://www.bls.gov/home.htm
Symbol : CES0000000001

## モデル構築

### 単回帰モデル

ADP雇用統計は実際の米国雇用統計の2日前に発表されます。
運用上この変数をモデルに入れてしまうと2日前にならないと予測値が出せません。
しかし、ADP雇用統計は相応に実際の雇用統計と相関が高く、先行指標として注目される指標でもあります。

まずは、これだけで単回帰モデルを作ってみます。
回帰モデルは単純で、$t$時点の雇用統計値を$t$時点のADPによって予測します。
(実際は時点がADPの方が早いですが、便宜上同時点としています。)

```math
LaborStat_t = ADP_t + e_t
```


一点注意すべき点は、$t-1$時点までの情報を使ってモデルを構築し、$t$時点の情報を使って予測を行うという点です。

```r
data2_lm ; #LaborStatとADPを格納したDataFrame
pred_vec <- rep(NA,12);

for(i in 12:(nrow(data2_lm)-1)){
  tmp_lm <- data2_lm[1:i,];
  lmret <- lm(LaborStat ~ ., tmp_lm);
  pred <- predict(lmret,data2_lm[i+1,]);
  pred_vec <- c(pred_vec,pred);
}
```

得てして面白みのないモデルですが、

![Rplot02_v0.jpeg](https://qiita-image-store.s3.amazonaws.com/0/54539/d549356d-48c2-deac-e71a-a6037f1d65cd.jpeg "Rplot02_v0.jpeg")

```r
print(cor(gdata$LaborStat,gdata$LaborStatPred)) #0.9459864
```

ん？おかしい。いけすぎてる。

```r
gdata2 <- gdata[gdata$date>=as.Date("2010-01-01"),];
print(cor(gdata2$LaborStat,gdata2$LaborStatPred)) #0.5981159
```

思った通り、リーマンショックの時のデータが相関を底上げをしてしまっているようです。

### 重回帰モデル

以下を使って重回帰モデルを構築してみます。

* ニューヨーク連銀製造業景気指数の「雇用」
* リッチモンド連銀製造業景況指数の「雇用」
* ＡＤＰ雇用統計
* Dow Jones月次変化率 (SPの代わり)
* 独自指標1個

回帰は基本同じですが、データが追加されたことに加えて推定方法が若干変わっています。

```r
library(dplyr)
quantile_data <- function(y,upper=0.9,lower=0.1){
  qt <- percent_rank(y);
  idx <- (qt < upper) & (qt > lower)
  which(idx);
}

data2_lm;
pred_vec <- rep(NA,12);

for(i in 12:(nrow(data2_lm)-1)){
  tmp_lm <- data2_lm[1:i,];
  idx <- quantile_data(tmp_lm$LaborStat,upper = 0.8,lower = 0.1);
  tmp_lm_idx <- tmp_lm[idx,];  
  lmret <- lm(LaborStat ~ .,tmp_lm_idx,weights = exp(seq(1,length.out = nrow(tmp_lm_idx),by=0.01)));
  slmret <- step(lmret,trace = F);
  pred <- predict(slmret,data2_lm[i+1,]);
  pred_vec <- c(pred_vec,pred);
}
```

モデル構築は全体のデータではなく、一部のデータを使用しています。
加えてstep wiseを行い、変数を逐次的に選択していきます。
また、回帰は直近を重くするように重みをつけ推定を行います。

![Rplot02_v10.jpeg](https://qiita-image-store.s3.amazonaws.com/0/54539/4e90948a-0a74-ccaf-f716-c5669c6bd7e6.jpeg "Rplot02_v10.jpeg")


```r
print(cor(gdata$LaborStat,gdata$LaborStatPred)) #0.9306076
```

```r
gdata2 <- gdata[gdata$date>=as.Date("2010-01-01"),];
print(cor(gdata2$LaborStat,gdata2$LaborStatPred)) #0.6119195
```

## とりあえずのまとめ

* リーマン時は結果の評価に使うと偽相関に陥るので、やめた方がいい。
* 変数を増やすと確かによくなるが、まだ不十分なレベル。
* そもそも雇用統計の結果が正規分布しているのか？
* 逆に変数てんこ盛りで機械学習のフレームワークに入れた方がいいかもしれない。(SVRとか)
* 完成までの道のりは長い。