要約
需要曲線を観測したデータから導くのは難しいです。なぜなら、観測できるのは、常に変化する需給曲線の交点だからです。
例えば、価格と量の関係性を分析しても、下図のように全く見当違いの結果を導きます。
この問題を解決するために、ライト親子によって考案されたのが操作変数法です。
本稿では、操作変数法が最小二乗法より優れていることを、シミュレーションで確かめました。
全体像
（1）問題の構造
問いとして「需要曲線を推定するシミュレーションをした際、操作変数法の方が最小二乗法より優れているか？」を設定します。この問いは以下の要素に分解できます。
①操作変数法と最小二乗法とは？
②シミュレーションの目的は？
③どうなれば成功か？
④どんな状況を想定するか？
⑤シミュレーションの結果は何か？
⑥それは成功か？
（2）前提の選択
①〜④を前提とします。
（3）論点の選択
⑤と⑥を論点します。
（4）付録一覧
長くなっちゃうので「小ネタ」と「Rコード」は付録に回しました。
前提
操作変数法
操作変数モデルについては「操作変数モデルとは？」の記事をご覧ください。
操作変数モデルをイメージ図で説明するとこんな感じです↓
重回帰モデルと比較するとわかりやすいかもしれません。↓
最小二乗法
本稿における最小二乗法とは、通常の単回帰分析を意味します。詳しくは「単回帰モデルとは?」「単回帰分析における最小二乗(OLS)推定量とは？」をご覧ください。
シミュレーションの目的
シミュレーションでは
・最小二乗法では需要曲線の傾きを上手く推定できない
・操作変数法では需要曲線の傾きを上手く推定できる
ということを示したいです。そのために、
・ある特定の1通りの場合で↑の結果になるのか
・1万通りを総合した場合で考えても↑の結果になるのか
を検証します。
シミュレーションの成功とは
「ある特定の1通りの場合」では、最小二乗推定値より、操作変数推定値の方が真の値に近ければ成功です。
「1万通りを総合した場合」では、最小二乗推定値より、操作変数推定値の分布の方が真の値に近ければ成功です。定量的には、1万通りの推定値の平均値で、推定の良し悪しを評価できます。
想定する状況
話を想像しやすくするために、具体的なストーリーを作ってみます。
基本モデル
マーシャル高校（※付録へ）では、おにぎりの価格と量がその日の需要曲線と供給曲線の交点に決まります。
生徒は安いとおにぎりをたくさん買いますし、食堂のおばちゃんは高いとたくさん作ろうとします。価格はワルラス校長が決めます。ワルラス校長（※付録へ）は、需要曲線と供給曲線を直感できる超能力者です。
基本は1個100円で、1日100個売れます。
需要ショック、供給ショックの導入
需給曲線は毎日同じではありません。毎日何かしらかの需要ショック、供給ショックがあります。
$$需要曲線：Q_{需要}=-P+200+Shock_{需要}$$
$$供給曲線：Q_{供給}=P+Shock_{供給}$$
$$需要Q_{需要}、供給Q_{供給}、価格P、ショックShock$$
需要ショックは、期待値0、標準偏差20の正規分布に従います。供給ショックは、期待値0、標準偏差15の正規分布に従います。このデータはRで生成できます。
#Rコードは付録：Rコード「シミュレーション・データの生成」へ供給ショックは、ある操作変数「Time」と相関しています。具体的な内容は、後述します。
ライト君
ライト君は以上の前提を知らず、100日分の観測データしか入手できません。
その中で授業レポート課題として「おにぎりの需要曲線を推定する」に取り組もうとしています。
方法
（１）ある特定の1通りの場合
ライト君はサンプル・サイズ100のデータから需要曲線を推定できます。一度だけシミュレーション・データを作って、推定できた結果が「ある特定の1通りの場合」です。
最小二乗推定値と操作変数推定値のどちらが、真の値に近いのか確かめます。
（２）1万通りを総合した場合
神の立場にいる我々は、シミュレーション・データと推定自体は何度も実行することができます。そこで、1万通りのデータ、1万通りのライト君、1万通りの推定結果を作ることができます。
最小二乗推定値と操作変数推定値の分布のどちらが、真の値に近いのか確かめます。定量的には、1万通りの推定値の平均値で、推定の良し悪しを評価します。
結果
（１）ライト君視点
ライト君はおにぎりの需要曲線を知りたいと考えました。
最小二乗法
「おにぎり価格とおにぎり量をプロットして、最小二乗法で単回帰分析すればいいじゃん！ Rでやってみよう！」
#Rコードは付録：Rコード「ライト君の単回帰分析」へ「お、結果が出たぞ。あれ、なんで価格1円が上がると、おにりぎの需要が0.6個増えるんだ？ しかも、1%有意だ。おかしいぞ」
$$(需要)=0.619(価格)+39.267+誤差項$$
「散布図を作ってみよう！ あれ、右上がりの需要曲線になってる？ 我が校のおにぎりはギッフェン財なの？ いやそんなはずはない！」
#Rコードは付録：Rコード「ライト君の散布図」へ
「あ、そうか。この点ってのは、需要曲線の上に現れるのではなくて、需要曲線と供給曲線の交点に現れるから、これを結んでも需要曲線にならないのか！」
操作変数法
ライト君は悩んでいると、食堂のおばちゃんが裏で怒られているのを見つけました。
上司「あなたね。毎日来る時間がバラバラすぎるのよ。早く来るとおにぎりたくさん作るし、遅く来るとあまり作らないのは良くないわ。もっとモチベーションを一定に保ってちょうだい。」
おばちゃん「すみません」
そのとき、ライト君は思いつきました。
「おばちゃんの出勤時間を操作変数にすればいいのでは！？ おばちゃんの出勤時間は、価格に与える関連性を持つけれども、生徒の需要とは無相関だから外生性がある。操作変数の2条件を満たしているぞ」
「目的変数を需要(Quantity)、説明変数を価格(Price)、操作変数を『おばちゃんが学校に来たのは始業何分前か(Time)』とすると
$$(Quantity)=\alpha(Price)+\beta+誤差項$$
$$操作変数：Time$$
だね！よし。早速、操作変数法で推定してみよう」
#Rコードは付録：Rコード「ライト君の操作変数法」へ「結果が出たぞ！ お！ ちゃんと価格が1円上がると需要が減る結果になってる。これで『右下がりの需要曲線になってるな』」
「ただ、10%有意かあ。サンプル・サイズがちょっと足りなかったかな… まあ、いいや結論としては
$$(需要)=-1.493(価格)+252.874+誤差項$$
です。最初に出した最小二乗法の結果
$$(需要)=0.619(価格)+39.267+誤差項$$
とはだいぶ違うね〜」
評価
答えは
$$(需要)=-1(価格)+200+誤差項$$
です。つまり、
$$真の値：-１$$
$$最小二乗推定値：0.619$$
$$操作変数推定値：-1.493$$
となります。1通りの場合、操作変数推定値に軍配が上がりました。
（２）1万人のライト君
次に、1万回シミュレーションデータを生成し、1万回推定した値についてみていきましょう。
#Rコードは付録：Rコード「1万回の操作変数法」へサンプル・サイズ100
サンプル・サイズ100で推定する場合を1万回試した場合、推定値のヒストグラムは次のようになりました。
これを見ると、最小二乗法によって推定すると、かなりバイアスがある一方で、操作変数法によって推定するとバイアスは比較的に小さいことがわかります。
操作変数推定値で得られる値の平均はー1.14でした。したがって、真の値ー１に近い値が得られたはずでした。とはいえ、得られうる推定値の最小値はー12でかなり誤差のある推定もあり得ました。
サンプル・サイズ300
実は、操作変数法には一致性はありますが、不偏性はありません。詳しくは「操作変数推定量に、不偏性はないが、一致性はある」をご覧ください。
そこで、サンプル・サイズを300に増やしてみます。
操作変数推定値で得られる値の平均はー1.03でした。推定値の最小値はー3.5で誤差はかなり縮小してます。
サンプル・サイズ1000
サンプル・サイズ1000では、操作変数推定値で得られる値の平均はー1.01でした。推定値の最小値はー1.65です。
評価
分布の形状では、サンプル・サイズが100、300、1000のときでも、操作変数推定値に軍配が上がりました。
推定値の平均値でも、その結果は変わりません。例えば、サンプル・サイズが100のとき、
$$真の値：-１$$
$$最小二乗推定値の平均：0.60$$
$$操作変数推定値の平均：ー1.14$$
でした。
考察
（1）結論
結論は「需要曲線を推定するシミュレーションをすると、操作変数法の方が最小二乗法より優れている」です。
（2）妥当性評価
前提評価
シミュレーションの目的、成功の判断基準、需給曲線について明確に示したのは、Goodです。
需要曲線・供給曲線の傾きが変わらないことを暗黙のうちに仮定していた点がBadです。
方法評価
1通りの場合、1万通りの場合について明確にシミュレーション方法を示していてGoodです。Rコードも付録についている点もGoodです。
結論評価
内生性がある場合、操作変数法が最小二乗法より優れているという結果が導けたのがGoodです。
かなり単純な場合を想定しているので、現実でも同じような結果が導けるかについては言及できないのがBadです。
（3）意義
操作変数法について理解が深まりました。
付録：小ネタ
マーシャル高校は、アルフレッド・マーシャル（イギリス、1842年-1924年）が元ネタです。”Principles of Economics”(1890年)で、消費理論と生産理論を統合しました。このとき、マーシャルは、需要曲線・供給曲線の図を用いました。実は需要曲線は限界効用理論という当時の消費理論、供給曲線は限界費用理論という当時の生産理論を意味しています。マーシャルは「価格を需要と供給のどちらが決めるのかは、ハサミの上の歯と下の歯のどちらが紙を切っているのかを問うのと同じくらい無意味だ」のような名言を残しています。「価格の説明には、需要と供給のどちらも不可欠だ」という意味ですね。
ワルラス校長は、ワルラスの競り人が元ネタです。需給曲線の交点で均衡します。ところで、少し経済学を勉強したことがある人は、プライステイカーの仮定を聞いたことがあると思います。これは消費者も生産者も、自分では価格を変えられず、価格を受け入れるしかない存在という仮定です。となると、市場の誰もが現状の価格を受け入れるしかなく、均衡していようが、不均衡だろうが、価格は変わらないわけです。おかしい！ そんな矛盾に応えるために作られたのがワルラスの競り人です。こいつは価格ごとの需要と供給を聞いて、両者が一致したときに交渉を成立させます。ん？ 「現実にそんなのはいないだろ」って？ それは言っちゃいけないお約束です。
ライト君は、フィリップ・ライト(父)とシューアル・ライト(子)の親子が元ネタです。1928年にフィリップ・ライトが発表した著書『動物性油と植物性油への関税』で操作変数法に書かれていました。もっというと、操作変数法について書かれていたのは付録Bでした。さて、付録Bは誰が書いたのでしょうか？ 父のフィリップ・ライトは研究者としてパッとせず、息子のシューアル・ライトは卓越した人口遺伝学者・統計学者だったので、「父親と息子のどちらが付録Bを書いたんだ？」問題が後に起こりました。どちらが付録Bを書いたのかは今となってはわかりません。
付録：Rコード
シミュレーションデータの生成
sample_size <- 100 #サンプルサイズの決定
Shock_Demand <- rnorm(n=sample_size, mean =0, sd=20) #需要ショック
Shock_Suply <- rnorm(n=sample_size, mean =0, sd=10) #供給ショック
Price <- 100 + 0.5*(Shock_Demand - Shock_Suply) #価格決定
Quantity <- 100 + 0.5*(Shock_Demand + Shock_Suply) #個数決定
U <- rnorm(n=sample_size, mean =-30, sd=5) #誤差項決定
Time <- Shock_Suply + U #おばちゃんの出勤時間
Data0 <- data.frame(Price,Quantity,Time,Shock_Demand,Shock_Suply,U) #データフレーム化ライト君の単回帰分析
OLS_reg <- lm(Quantity ~ Price ,data=Data0) #単回帰分析
summary(OLS_reg) #詳細閲覧
library(stargazer) #回帰分析表のためのパッケージ呼び出し
stargazer(OLS_reg,
title = "おにぎりの需要曲線推定(メモ：失敗した！)", #図表タイトル
covariate.labels =c("おにぎり価格","定数項"), #説明変数の名前付
dep.var.caption = "目的変数",
dep.var.labels = "おにぎり個数", #目的変数の名前付
type="text", #テキスト・タイプでRStudioにて出力
omit.stat = c("adj.rsq","ser","f"))ライト君の散布図
plot(Quantity,Price) #散布図
library(ggplot2) #綺麗な散布図のためのパッケージ呼び出し
ggplot(Data0, aes(x = Quantity,y=Price))+ #データ指定
geom_point()+ #散布図
scale_x_continuous(limits = c(70, 140),breaks = seq(70,140,10))+ #縦軸の目盛り調整
scale_y_continuous(limits = c(70, 140),breaks = seq(70,140,10))+ #横軸の目盛り調整
stat_smooth(method = "lm")+ #回帰直線の描画
xlab("おにぎり個数") + #横軸ラベル
ylab("おにぎり価格") + #縦軸ラベル
theme_grey(base_family = "HiraKakuPro-W3") #日本語指定ライト君の操作変数法
#操作変数法のためのパッケージ呼び出し
library(AER)
#操作変数法の実行
IV_reg <- ivreg(data = Data0, Quantity ~ Price | Time)
#ivreg・・・操作変数法
#data = Data0はData0というデータフレームを使ったということ
#Y ~ X | Z はYが目的変数、Xが説明変数、Zが操作変数
#結果の出力
summary(IV_reg)
#綺麗な回帰表の出力
stargazer(IV_reg,
title = "おにぎりの需要曲線推定(メモ：操作変数法)",
covariate.labels =c("おにぎり価格","定数項"),
dep.var.caption = "目的変数",
dep.var.labels = "おにぎり個数",
type="text",
notes="操作変数としておばちゃんの出勤時間を利用 ",
omit.stat = c("rsq","adj.rsq","ser","f"))1万回の操作変数法
library(AER) #操作変数法パッケージ呼び出し
library(ggplot2) #作図パッケージ呼び出し
sample_size <- 100 #サンプル・サイズ
trial_number <- 10000 #試行回数=1万
OLS <- rep(NA, trial_number)
IV <- rep(NA, trial_number)
for(j in 1:trial_number){
Shock_Demand <- rnorm(n=sample_size, mean =0, sd=20)
Shock_Suply <- rnorm(n=sample_size, mean =0, sd=10)
Price <- 100 + 0.5*(Shock_Demand - Shock_Suply)
Quantity <- 100 + 0.5*(Shock_Demand + Shock_Suply)
U <- rnorm(n=sample_size, mean =-30, sd=5)
Time <- Shock_Suply + U
data0<- data.frame(Price,Quantity,Time,Shock_Demand,Shock_Suply,U)
OLS_reg <- lm(Quantity ~ Price ,data=data0)
IV_reg <- ivreg(data = data0,Quantity ~ Price | Time)
OLS[j] <- coefficients(summary(OLS_reg))[2,1]
IV[j] <- coefficients(summary(IV_reg))[2,1]
}
data_OLS_IV <- data.frame(OLS,IV)
ggplot(data_OLS_IV) +
geom_histogram(aes(x = IV),colour="skyblue3", fill = "skyblue2",alpha=1,bins=150)+
geom_histogram(aes(x = OLS),colour="pink3", fill = "pink1",alpha=1,bins=150)+
geom_vline(xintercept =-1 )+
geom_vline(xintercept =mean(IV), colour="blue")+
geom_vline(xintercept =mean(OLS), colour="red")+
xlab("操作変数による推定値(青)、最小二乗法による推定値(赤)、パラメーター(黒)") +
ylab ("度数" )+
theme_grey(base_family = "HiraKakuPro-W3")
カテゴリー