【統計】【R】IPW推定量を導出してみる その３

前回の記事（その１、その２）で紹介したIPW推定量のRコードを紹介します．

次のような２つのステップを踏んで、IPW推定量を計算します．
傾向スコアを推定
傾向スコアの推定には、ロジスティクス回帰分析を用います．割り当て変数を従属変数、共変量を説明変数にして分析を行います．
IPW推定量を計算
前回の記事より、IPW推定量を計算します．

サンプルデータはRのlalondeというデータを用います．
このデータは、職業支援プログラムの賃金（78年の賃金：re78）への効果をはかるもので、訓練を受けるグループ（処置群：treat=1）と受けないグループ（対照群：treat=0）間の賃金差を調べています．

	#estimation for propensity score
	data(lalonde)
	logi <- glm(treat ~ age + educ + black + hisp + married + nodegr +
	re74 + re75 + u74 + u75, family = binomial, data = lalonde)
	#IPW estimation
	IPW1 <- sum( (lalonde$re78/logi$fitted)[lalonde$treat == 1]) /
	sum(1/logi$fitted[lalonde$treat == 1])
	IPW0 <- sum( (lalonde$re78/(1-logi$fitted))[lalonde$treat == 0]) /
	sum(1/(1-logi$fitted)[lalonde$treat == 0])

view raw IPW estimator.R hosted with ❤ by GitHub

応援よろしくです．

【統計】IPW推定量を導出してみる その２

前回の記事の続きです．Horvitz-Thompson型推定量の推定値を求めます．


Horvitz-Thompson型推定量
「強く無視できる割当条件」\[ (y_{0}, y_{1}) \perp z | \mathbf{x} \]が成立するとき、
\[ E\left[\frac{y_{1}z}{e(\mathbf{x})}\right] = E\left[y_{1}\right] .\]

期待値の推定量は平均値をとることが一般的なので、


\begin{equation*}
\begin{split}

E \left[ \frac{y_{1}z}{e} \right]
&= \frac{1}{n} \sum_{i=1}^{N} \frac{z_{i}y_{i1}}{e_{i}} \\
&= \frac{\sum_{i=1}^{N} \frac{z_{i}y_{i1}}{e_{i}}}{\sum_{i=1}^{N} \frac{z_{i}}{e_{i}}}
\end{split}

\end{equation*}

1行目から2行目へは、次のような関係を用いています．


\begin{equation*}
\begin{split}

E \left[ \frac{z}{e} \right]
&= E_{\mathbf{x}} \left[ E \left[ \frac{z}{e} | \mathbf{x} \right] \right] \\
&= E_{\mathbf{x}} \left[ \frac{1}{e}E \left[ z | \mathbf{x} \right] \right] = 1．

\end{split}
\end{equation*}

ところで、
\[ E \left[ \frac{z}{e} \right] \approx \frac{1}{n} \sum_{i=1}^{N} \frac{z_{i}}{e_{i}}\]
であるので、
\[ \frac{1}{n} \sum_{i=1}^{N} \frac{z_{i}}{e_{i}} \approx 1 \leftrightarrow
\sum_{i=1}^{N} \frac{z_{i}}{e_{i}} \approx n． \]

次回は、この推定量のRコードを考えてみます．

参考文献

応援よろしくです．

【統計】IPW推定量を導出してみる その１

前回の続きで、今回はRubin(1985)の提案している「傾向スコアによる重み付け推定法」を導出しようと思います．

• $z$：割付け変数
• $x$：共変量
• $y$：従属変数

IPW推定量は、傾向スコアの逆数の重み付け平均でHorvitz-Thompson型推定量と呼ばれています．この推定量は通常の単純平均とは異なり、共変量を調整した平均です．この推定量のいいところは、平均処遇効果と周辺平均を計算することができることです．

マッチングや重回帰分析による共変量調整は、平均処遇効果を計算することはできますが、割り付け変数ごとに周辺平均を計算することはできません．しかし、この推定量を用いれば、周辺平均を計算することができます．

Horvitz-Thompson型推定量
「強く無視できる割当条件」\[ (y_{0}, y_{1}) \perp z | \mathbf{x} \]が成立するとき、
\[ E\left[\frac{y_{1}z}{e(\mathbf{x})}\right] = E\left[y_{1}\right] .\]

証明
まず傾向スコアの定義から
\[ e_{i} \equiv Pr(z_{i} = 1|\mathbf{x_{i}})  = E(z_{i}|\mathbf{x_{i}}) \]
と書くことができます．以下では対象者番号のｉを省略しています．

\begin{equation*}
\begin{split}

E \left[ \frac{y_{1}z}{e} \right]
&= E_{ \mathbf{x}} \left[ E\left( \frac{y_{1}z}{e}| \mathbf{x} \right) \right] \\
&= E_{ \mathbf{x}} \left[ E\left(y_{1}|\mathbf{x}\right) E\left( \frac{z}{e}|\mathbf{x}\right) \right] \\
&= E_{ \mathbf{x}} \left[ \frac{1}{e} E\left( y_{1}|\mathbf{x} \right) E\left( z | \mathbf{x} \right) \right] = E\left( y_{1} \right)
\end{split}

\end{equation*}


1行目から2行目へは割付が強い意味で無視できるという条件を用いています（$\mathbf{x}$が与えられた下では、$y$と$z$は条件付き独立）Q.E.D.

上の方法と同様にして、$E[y_{0}]$を求めることができます．


\begin{equation*}
\begin{split}

E \left[ \frac{y_{0}(1-z)}{1-e} \right]
&= E_{ \mathbf{x}} \left[ E\left( \frac{y_{0}z}{1-e}| \mathbf{x} \right) \right] \\
&= E_{ \mathbf{x}} \left[ E\left(y_{0}|\mathbf{x}\right) E\left( \frac{1-z}{1-e}|\mathbf{x}\right) \right] \\
&= E_{ \mathbf{x}} \left[ \frac{1}{1-e} E\left( y_{0}|\mathbf{x} \right) E\left( 1-z | \mathbf{x} \right) \right] = E\left( y_{0} \right)
\end{split}

\end{equation*}



つまり、真の傾向スコアがわかっていて、強く無視できる割当ての条件が成立するならば、対照群・処置群それぞれの周辺期待値を求めるには、傾向スコアによって重み付けされた期待値を求めればよい、ということになります．

参考文献

【統計】傾向スコア

傾向スコアが因果推論でよく使われているということで勉強中…

ということで、「傾向スコアを用いた共変量調整による因果効果の推定と臨床医学・薬学・公衆衛生分野での応用について」の前半部分をまとめました．

共変量の影響を除去しなければならない！

たとえば、ある薬を飲むこと（従属変数）で病気が早く治る（独立変数）か知りたいとき、投薬の有無を操作することができないために、薬を飲んだ人と飲まなかった人に分けて治癒日数の平均の差を調べることになります．しかし、この方法では治癒日数に影響を与えそうな病気の症状や年齢などの共変量の影響が除去されないために、ある薬の病気を早く治す因果効果を知ることができなくなってしまう！
じゃどうすればいいのか？？ということで、共変量の調整方法として以下の３つが考えられてきました．

既存の共変量調整

 既存の共変量調整法は３つあります．

1. 均衡化
共変量が同じになるペアを作ること
例えば、病気の症状を重いと軽いに分けること

2. 恒常化・限定
共変量がある値で一緒の被験者のみに限定して分析を行うこと

3. 統計的調整
1. マッチング
2. 層別解析
3. 共分散分析

最後の３．共分散分析は、共変量が独立変数と従属変数に影響を与えてしまうなら、それをモデルに含めて考えればいいという方法になります．

具体例
Cox回帰、パス解析、グラフィカルモデリング、構造方程式モデリング…
問題点
しかし、この方法でもやはり問題点があります。
線形や二次関数など、従属変数と共変量の関係をモデル化しないといけないということです．もし真のモデルが二次関数のときに線形を仮定してしまうと誤った結論が得られてしまいます．また、回帰係数自体が、因果効果と等しくはならないということです．

傾向スコアの利用

そこで、こうした問題点を克服するために、”セミパラメトリックな”共変量調整法として傾向スコアを用いた解析法が注目されています．
定義
今、薬を飲んだ人をz＝１、薬を飲まなかった人をz＝０のようにzを割付変数とします．また共変量ベクトルをxとすると、傾向スコアは群１（薬を飲んだ人）へ割付けられる確率として定義されます．添字のiは第i被験者を表しています．

\[ e_{i} = Pr(z_{i} = 1|\mathbf{x_{i}}) \]

使用のための条件
ただし、この傾向スコアを使うためには条件が必要です．それが、割当てが「強い意味で無視可能」であるという条件です．この条件は次の３つに分けられます．

1. 観測された共変量によって割付が説明できるということ
割付は共変量のみから影響を受けるということ
たとえば、薬を飲むかどうかに割付けられるかは、共変量である病気の症状によって説明できるということです．


2. 観測されていない共変量が割付に影響を与えないということ
観測したもの以外に共変量が存在しないということです．
傾向スコアを用いた解析では、アンケート調査で共変量をもれなく聞いている必要があります．


3. 従属変数の値によって、割付が決定されないということ
共変量を与えた下では、病気の治癒日数である従属変数が投薬をするかしないかの割当てには影響を与えないということです．つまり、以下のような条件付き独立が成立していることが必要です．

傾向スコアの具体的な解析方法
以上のような条件が成立するときに、傾向スコアを用いた解析を行うことができます．
その解析方法は二段階に分かれています．

1. 傾向スコアの推定
ロジスティック回帰かノンパラメトリック回帰で推定されます．


2. 傾向スコアを用いた調整
1. マッチング
2. 層別
3. 共分散分析
4. 重み付き平均
マッチング・層別・共分散分析以後、重み付け平均の方法が提案されていて、次第に利用されてきています．

次回は、重み付け平均を考えてみます．

【R】【ファイナンス】ヒストリカルボラティリティを計算してみる

ヒストリカルボラティリティ（HV）を計算できるプログラムをRで書いていました．

20日間で計測することが多いようなので、計算式は以下のようにしています．

HV = 20日間の対数収益率の標準偏差 × 250の平方根

HVは過去の株価の変動が大きくなれば大きくなり、一定の割合ならば小さくなります．

HV導出に用いたデータはyahoo! financeから取得しています．
データ：Nikkei 225

サンプル期間：2012年1月4日〜5月14日

Rスクリプト

	#20日間のHVをまとめて表示
	#yahoo! Financeから情報取得
	test2 <- read.csv("table.csv", header = T)
	attach(test2)
	#標本分散
	variance <- function(xx){
	(var(xx, na.rm = T) * length(xx - 1)) / length(xx)
	}
	#対数収益率を計算
	rat <- c()
	for(ii in 2:length(Date)){
	rat[ii] <- log(Close[ii]/Close[ii - 1]) * 100
	}
	#HVを20日間計算
	dd <- c()
	for(ii in 1:length(Date)){
	diff <- length(Date) - 19
	cc <- c()
	cc <- rat[(ii+1):(ii+20)]
	if(ii <= diff){
	dd[ii] <- sqrt(variance(cc)) * sqrt(250)
	}
	else {
	break
	}
	}
	#作成した対数収益率とHVをデータフレームに追加
	temp <- rev(dd)
	temp <- c(temp, rep(NA, 19))
	temp <- rev(temp)
	test2 <- transform(test2, "Return"=rat ,"HV"=temp)
	#対数収益率とHVをプロット
	par(mar=c(5, 5, 4, 5) + 0.1)
	plot(test2$Return, type="l", ylim=c(-4,4), xlab="", ylab="", col="red")
	par(new=T)
	plot(test2$HV, type="l", axes=F, xlab="", ylab="", col="blue")
	axis(4)
	mtext("HV", side=4, line=3)
	title("Rate of Log Return and HV", xlab="Time", ylab="Return")

view raw historical volatility.R hosted with ❤ by GitHub

図

赤線が対数収益率、青線がHVを表しています．4月の半ば頃からHVが高くなっている様子が伺えます．

参考書


「R言語逆引きハンドブック」はRコードを書くときにお世話になりました．

応援よろしくです．

【論文】【ファイナンス】【言語処理】ツイッターで株式市場予測？！

ツイッターから感情を抽出し、ダウ工業平均の変動を予測した論文Twitter mood predicts the stock market.
ツイートから社会全体の雰囲気を表す感情指標を作成し、それを使って株式市場を予測しようという論文．

---

【概要】

• Twitterの投稿の語彙から抽出された社会全体の感情を表す指標が、DJIAの変動を87.6%の精度で予測した．
• 3〜4日前のCalmという感情の変動がDJIAの変動とほぼ一致する．

【解析の流れ】

• OFとGPOMSを用いて、7つの感情に分類する．Yahoo! FinanceからDJIAを取得する．
• OFとGPOMSを比較できるように平均0分散1のZ値に標準化する．
\[ \mathbf{Z_{X_{t}}} = \frac{X_{t} - \bar{x}(X_{t \pm k})}{\sigma(X_{t \pm k})} \]
• グレンジャー因果性テストでOF¥&GPOMSがDJIAの予測に使えるかを検証する．
• SOFNNを用いて非線形のDJIA予測モデルの予測精度が感情因子の組み合わせによって向上するかどうかを確認する．

【OFとGPOMS】

• OF
• センチメント分析のフリーソフトで、文章の感情を特定できる．
• 2718個のpositive語と4912個のnegetive語からなるOF辞書を使って、ある日のツイートの感情をpositiveとnegativeに分類し、それぞれの感情が投稿されるとそれぞれに1だけ増やし、その比率を計算する．
• p/nの2つだけであるので、人間の複雑な感情を把握できているとは言い難い．


• GPOMS
• ツイートのテキストデータをCalm, Alert, Sure, Vital, Kind, Happy に分類する．
• 2006年にGoogleがwebサイトから抽出した1兆語もの単語を利用し、25億に及ぶ4-/5-gramsの解析から作成された．

【グレンジャー因果性】

• 3日前のCalmの変動が現在のDJIA変動と同じような動きをしていることから、Calmの時系列データがほかの感情データと比べて予測の情報を含んでいると言える．
• ただし、10月13日の連邦準備制度理事会のアナウンスの結果、Calmのデータからは予測できないDJIAの変動が生じている．

【SOFNN】

• OFのp/nは予測精度を向上させることはない．
• Calmがもっとも高い予測精度（株価の上下変動の予測）を記録した（86.7¥%）．
• CalmとSureやCalmとVitalを組み合わせたものは、将来のDJIAの情報を含まない．CalmとHappyの組み合わせは、MAPEをもっとも下げ、株価の上限変動の予測に関していい結果を残した（80.0¥%）．Happyはグレンジャー因果性のテストでは、DJIAの予測には役立たないという結果であったが、Calmと組み合わせるとDJIAの予測に役立つ．Fテストをすると、CalmとHappyを線形結合するとCalmだけの場合に比べて結果がよくないので、異なった感情因子には非線形の関係があることが確認される．

---

【感想】

• グレンジャー因果性を経済学の文脈で用いるとき、次数は1にすることが多いが、本論文では、7までテストしている．
• 収集したツイートは必ずしも市場の情報を含んでいるものとは言えない．社会全体の感情や気分が株式市場を予測できるかどうかが本論文の目的であるが、マーケットに関連する、あるいは投資家のクラスタの中で同様の分析をすると予測精度がさらに上昇しそう．
• FRBがBail-outのアナウンスをしたときには、感情データは株式市場を予測できていない．このあたりは経済学的な考え方で株価上昇の説明ができる．
• 感情データと株式市場には相関がありそうだが、因果関係があるかは分からない．
• パフォーマンスは？

応援よろしくです．

【お知らせ】【統計】勉強会をやろうと思っています

回帰分析だけでなく、統計分析を幅広く勉強しようと思っていたときに、友人に勧められて購入した本．



　現象のモデル化ということで回帰モデルはよく使いますが、識別／判別／分類はあまり使わないので、勉強になります．
　経済学では、特に回帰分析が多いな、というのが私のイメージです．ビジネスの分野では、識別／判別／分類などの手法がよく使われているそうです．

目次

1. はじめに
2. 線形回帰モデル
3. 非線形回帰モデル
4. ロジスティック回帰モデル
5. モデル評価基準
6. 判別分析
7. ベイズ判別
8. サポートベクターマシーン
9. 主成分分析
10. クラスター分析

付録ブートストラップ法
付録ラグランジュの未定乗数法
付録EMアルゴリズム

【お知らせ】

この本を使って、勉強会をしようと考えています．

本ブログを通じて、詳細をお知らせできればと思います．

場所は大阪です．

【気になる記事】【統計】WEBで読める統計関係の良質な資料 | など

最近の気になった記事を紹介しています．

WEBで読める統計関係の良質な資料

独学で勉強できる環境になっています．

統計関係で、この他にもお世話になっているサイトをこのブログでも紹介できればと思います．

グーグル、検索システムを改良へ　市場シェアの維持目指し過去最大の変更か

検索サイトからSMSに人が流れている中で、グーグルはその地位を守ることができるかでしょうか？

FC2ブログランキングに登録しました．

ポチッと押して応援よろしくお願いします．

【気になる動画】【動物】犬のような猫

猫も犬も大好きなので、犬的要素を兼ね備えた猫はサイコーです．



リンク先

Cat Does Dog Impression

【気になる記事】【ライフハック】これはマネしたい！模様替えの参考になりそうなパソコン周りの画像を集めました(*´ω｀*)｜など

気になった記事の紹介です．

これはマネしたい！模様替えの参考になりそうなパソコン周りの画像を集めました(*´ω｀*)

パソコン周りの模様替えにおすすめの記事です．
Macがかっこいいー！

旦那さんのやみつき油そば

自分で料理する訳ではないけど、クックパッドにはまっています．

【読書メモ】【行動経済学】【効用理論】誘惑される意志２

今回は、前回「誘惑される意志」の続きを書いていきたいと思います．




第二章　意思決定の科学の根底にある二律背反

本章では、中毒者の例を挙げながら効用理論の限界が記述されています．
今回の投稿ではそのことをメインで書いていきたい思います．

効用理論

• 効用理論のもとで、人は、選択肢それぞれに報酬があって、それを最大化しようとしている．
• 効用理論において不合理な行動は、主要な確率計算の失敗が原因．（非合理性＝期待収入の失敗）
• 効用理論に説明できること：低金利の預金に手をつけたくなるからといって、高金利のローンに手をだすような行動がどのくらい余計なコストをもたらすか
• 説明できないこと：薬物中毒、悪習、ギャンブルのような変な行動をなぜやめないのか？

---

【感想】

• 本書にも書いているように「合理的な中毒者」に対する批判
• 今まで、何も気にせず効用最大化問題を解いてきたけど、こういう理解も必要なんだと反省
• 勉強の合間に悪習と思いながらついついツイッターを開いてしまう、行動をどうにかしたい笑”

【読書メモ】【行動経済学】誘惑される意志

去年の夏に図書館で借りて、去年の秋に購入した「誘惑される意志」．ゆっくり読んでもあまり理解できずに本棚にお眠りになっていたけど、双曲線割引についてもっと理解せねばということで、また読んでます．今度こそ！という思いで、ブログに書いて何とか続けようと努力しているわけで…

章ごとに感想とかまとめとかを書いていきます．

筆者は精神科医のジョージ・エインズリー．

行動経済学の論文にはよく登場するみたいで．
秋に読んだBenabou and Tirole(2004)にもばっちり引用されてました．
始めて勉強したときは、こんなことを経済学で説明できるのか、と思いました．

第一章　はじめに

副題にもあるように、人はなぜ自滅的行動をするのか、っていうのが本書のテーマ．
第一章には自滅的な行動をしてしまう、ことについての観察や研究が書かれています．
また、その学問領域も紹介されています．

自滅的な行動の小史／自滅的な行動の研究方法
経済学専攻なので、主に経済学や効用理論のところを重点的に見ているのですが、認知主義や精神分析、神経生物学などさまざまな切り口から人間の意志についての研究の概要が簡単にまとめられていて、とても興味深いです．

気になった経済学の説明は…
エインズリーの経済学の合理的意志決定の説明
合理的な意思決定：選択肢の熟考の結果／外部の財に基盤／新しい情報がない限り安定
非合理的な行動の理由：不適切な情報／将来の大幅な割引

「外部の財に基づかない（頭の中の内的な報酬と結びつかない）」と書いているけど、これってどういう意味でしょうか？

人間の行動は何らかの欲求や要求を満たすために起こっていて、その欲求（要求）の対象が「外部の報酬」にあれば外発的動機づけと呼び、欲求（要求）の対象が『内部の価値観（興味関心）』にあれば内発的動機づけと呼ぶそうです．

外的な、と言われると例えば政府の政策のようなものを想像します．
こんな行動をしたら報酬がもらえますよ、とか罰せられますよとか．
子ども手当なら（実際はそうじゃないかもしれないけど）子どもの数が増えればお金がたくさんもらえるから、子どもたくさん作りましょう！っていうのが、外的な報酬の例でしょうか．

一方、内的な、と言われると、よく学校教育で叫ばれているような「自律」とか「能動的に学びましょう」っていう言葉を連想します．
その行動自体に興味関心を抱き、目標を達成していくことに喜びを感じる、というのが内的な報酬のイメージになるでしょうか．

本書では、「自滅的な行動」の理由として、経済理論の不合理性の理由では不適切であることを紹介しています．

【お知らせ】【大学院】大学院受験

今日は、後輩の書いている大学院受験ブログを紹介します．

http://bei0135.blog79.fc2.com/

専攻は経済学で、京大／阪大／神大を受験してすべて合格したそうです．
（すごい！）

いろんな本を紹介しているので、
これから大学院（経済）を受験しようとしている人はぜひ、参考にしてみてください．

彼は、大学院の受験が終わったので、
この先、ブログの方向性がどうなるのか、結構楽しみにしてます．

【読書メモ】【統計】【R】R言語逆引きハンドブック

前回に引き続きRのことを書いています．


Rの書き方や関数はほとんどネットに落ちているので、困ることはほとんどありませんが、
辞書的なものも欲しいということで、かなり高かったのですが、R言語逆引きハンドブックを購入。

パソコンの隣に置いてけば安心です．

この本はOsaka.Rという統計とRの勉強会で紹介されていた本でもあります．

最近、経済学から離れて統計とプログラムがメインになっています．