第7回 パネルデータ解析演習

---

詳細はスライド http://user.keio.ac.jp/~katayama/lecture/panel.pdf を参照してください．

7.1 パネルデータにおける因果効果の推定

必要なパッケージの読み込み．必要に応じてインストールしてください．

library(dplyr)
library(tidyr)
library(plm)
library(wooldridge)

パネルデータ

• ユニットの特徴を一定期間に渡って繰り返し観測したデータ
• 次の例で見るようなロング型で表される．

パッケージwooldridge内のデータwagepanを読み込む．

• nr : ユニット識別番号
• lwage : 対数変換された賃金
• union : 労働組合に参加しているかどうか
• exper : 労働市場加入年数

data(wagepan)
wage1 <- wagepan[,c("nr","year","lwage","union","exper")]
head(wage1,10)

A data.frame: 10 × 5

	nr	year	lwage	union	exper
	<int>	<int>	<dbl>	<int>	<int>
1	13	1980	1.1975402	0	1
2	13	1981	1.8530600	1	2
3	13	1982	1.3444617	0	3
4	13	1983	1.4332134	0	4
5	13	1984	1.5681251	0	5
6	13	1985	1.6998910	0	6
7	13	1986	-0.7202626	0	7
8	13	1987	1.6691879	0	8
9	17	1980	1.6759624	0	4
10	17	1981	1.5183982	0	5

一般的なモデル $$ Y_{it}= \delta^{T}D_{it} + \theta^{T}X_{i} + u_{i} + \varepsilon_{it},\quad i=1,\dots, N; t=1,\dots T$$

• $N$ : ユニット数, $T$ : 時点数
• $Y_{it}\in\mathbb{R}$ : Outcome(1次元)
• $D_{it} \in \mathbb{R}^{K}$ : 興味のある変数($K$次元)
• $u_{i}$ : 未観測の交絡変数で，時点に依存しないunit-specificな変数
• $X_{i}$ : 観測される交絡変数で，unit-specific

Goal : $D$から$Y$への因果効果$\require{color}\textcolor{blue}{\delta}$を推定したい

• 簡単のため$X$は省略し，未観測の交絡変数$u$の存在は認める
• 観測データは各ユニット$i\in\{1,\dots, N\}$に対して$T$時点までのデータ

$$ \{(Y_{it}, D_{it}) : t=1,\dots, T\},\quad D_{it} = (D_{it1},\dots, D_{itK})^{T} $$

• 各ユニットはiid : $\{(Y_{i}, D_{i}, u_{i})\}_{i=1}^{n}$は$(Y, D, u)$の分布からのiidコピー．ただし，

$$ Y_{i}=\begin{pmatrix} Y_{i1}\\ \vdots \\ Y_{iT} \end{pmatrix},\;\; D_{i}=\begin{pmatrix} D_{i11} & \cdots & D_{i1K}\\ \vdots & \ddots & \vdots \\ D_{iT1} & \cdots & D_{iTK} \end{pmatrix} $$

Notes

• 未観測の交絡$u$が存在するケース $\rightarrow$ unobserved heterogeneity
  • $u$は固定効果(fixed-effects)とも呼ばれる
• 因果効果$\delta$はFixed-effects estimatorで推定可能
  • Within estimatorとも呼ばれる

7.1.1. Fixed-effects (within) estimator

考えるモデル $$ \begin{gather*} Y_{it} = \delta^{T}D_{it} + u_{i} + \varepsilon_{it} \\ \mbox{where}\quad E(\varepsilon_{i}\,|\,D_{i}, u_{i}) = 0_{T} \end{gather*} $$

Pooled OLS $$ \hat{\delta} = {\rm arg}\min_{b}\sum_{i=1}^{N}\sum_{t=1}^{T}(Y_{it} - b^{T}D_{it})^{2} $$

• 未観測の交絡変数を無視してOLSを適用したもの
• 交絡変数$u_{i}$を誤差項にpoolしている
• 興味のある変数$D_{i}$が$u_{i}$と相関 $\rightarrow$ バイアスが生じる(why?)

Fixed-effects (within) estimator $$ \require{color} (\hat{\delta}, \hat{u}_{1},\dots, \hat{u}_{N}) = \mathop{\rm argmin}_{b,m_1,\dots, m_N}\sum_{i=1}^{N}\sum_{t=1}^{T}(Y_{it} - b^{T}D_{it} - \textcolor{blue}{m_{i}})^{2} $$

• 交絡変数も最適化変数に入れて推定

First-order conditionsから$\delta$の推定量が陽に得られる :

$$ \begin{gather*} \hat{\delta} = \bigg(\sum_{i=1}^{N}\sum_{t=1}^{T}\ddot{D}_{it}\ddot{D}_{it}^{T}\bigg)^{-1} \sum_{i=1}^{N}\sum_{t=1}^{T} \ddot{D}_{it}\ddot{Y}_{it} \\ \mbox{where}\quad \ddot{D}_{it} = D_{it} - \bar{D}_{i},\quad \ddot{Y}_{it} = Y_{it} - \bar{Y}_{i},\quad \bar{Y}_{i} = \frac{1}{T}\sum_{t=1}^{T}Y_{it},\quad \bar{D}_{i} = \frac{1}{T}\sum_{t=1}^{T}D_{it} \end{gather*} $$

7.1.2 具体的なFixed-effects estimatorの計算

最初のwagepanのデータに対して，以下のモデルを想定する :

$$ {\tt lwage}_{it} = \delta_{1}{\tt union}_{it} + \delta_{2}{\tt exper}_{it} + u_{i} + \varepsilon_{it} $$

すなわち，$Y_{it} = {\tt lwage}_{it}$, $D_{it1} = {\tt union}_{it}, D_{it2} = {\tt exper}_{it}$と考えている． まずは$\ddot{Y}_{it} = Y_{it} - \bar{Y}_{i}$および$\ddot{D}_{it}=D_{it} - \bar{D}_{i}$をR上で計算する．上記のようなロングデータから$\bar{Y}_{i} = T^{-1}\sum_{t=1}^{T}Y_{it}$などの時点平均を計算するためには，各ユニットでグループ分けした上で，時点を渡って平均を取り，それを元の$Y_{it}$から除外する必要がある．

wage2 <- group_by(wage1,nr)  ##パッケージdplyr内のgroup_by関数を使う．変数nrによるグループ分け

見た目にはわからないが，上記のwage2はnrに基づいてグループ化されている．summarize関数と組み合わせると各グループにおける--この場合はユニットごとの--平均や標準偏差が計算できる．

head(summarize(wage2,mean=mean(union),sd=sd(union)))

A tibble: 6 × 3

nr	mean	sd
<int>	<dbl>	<dbl>
13	0.125	0.3535534
17	0.000	0.0000000
18	0.000	0.0000000
45	0.250	0.4629100
110	0.125	0.3535534
120	0.000	0.0000000

mutate関数と組み合わせて，$\ddot{Y}_{it}$および$\ddot{D}_{it1}, \ddot{D}_{it2}$の列を新たに追加する．グループ化済みのwage2に関する処理なので，mean関数はユニットごとに計算されていることに注意．

wage3 <- mutate(wage2,dlwage=lwage-mean(lwage),
                                dunion=union-mean(union), dexper=exper-mean(exper))
head(wage3,10)

A grouped_df: 10 × 8

nr	year	lwage	union	exper	dlwage	dunion	dexper
<int>	<int>	<dbl>	<int>	<int>	<dbl>	<dbl>	<dbl>
13	1980	1.1975402	0	1	-0.05811191	-0.125	-3.5
13	1981	1.8530600	1	2	0.59740793	0.875	-2.5
13	1982	1.3444617	0	3	0.08880960	-0.125	-1.5
13	1983	1.4332134	0	4	0.17756128	-0.125	-0.5
13	1984	1.5681251	0	5	0.31247305	-0.125	0.5
13	1985	1.6998910	0	6	0.44423889	-0.125	1.5
13	1986	-0.7202626	0	7	-1.97591466	-0.125	2.5
13	1987	1.6691879	0	8	0.41353583	-0.125	3.5
17	1980	1.6759624	0	4	0.03817607	0.000	-3.5
17	1981	1.5183982	0	5	-0.11938821	0.000	-2.5

全く同様のことは，以下のようにパイプ記号%>%を使ってもできる．ただしパッケージdplyrは必要である．次のコードは，wage3というオブジェクトに，wage1をgroup_by関数にかけて，その後にmutate関数で処理したものを代入せよというものである．ドットの部分は省略可能だが，その部分にパイプ記号の一つ前のオブジェクトが代入されると思えばok.

wage3 <- wage1 %>% group_by(.,nr) %>% 
  mutate(., dlwage=lwage-mean(lwage),dunion=union-mean(union),
                   dexper=exper-mean(exper))

これでFixed-effects estimatorを計算する準備ができた．その定義を思い出すと，

$$ \hat{\delta} = \bigg(\sum_{i=1}^{N}\sum_{t=1}^{T}\ddot{D}_{it}\ddot{D}_{it}^{T}\bigg)^{-1} \sum_{i=1}^{N}\sum_{t=1}^{T} \ddot{D}_{it}\ddot{Y}_{it} $$

である．ただし，$\ddot{D}_{it}$は一般に$K$次元ベクトルであることに注意．今回のケースでは$K=2$である．また，$nr=13$を$i=1$に，$year=1980$を$t=1$にそれぞれ対応させれば，

$$ \ddot{D}_{11} = \begin{pmatrix} -0.125 \\ -3.5 \end{pmatrix},\quad \ddot{Y}_{11}=-0.058 $$

である．以上を踏まえて$\hat{\delta}$を計算すると次のようになる． 適切にデータを作っておけば，lm関数でも計算できる．(why?)

##行列形式に変換．これをしておかないと，apply関数が適用できない場合がある(wage3はtibble形式のため)．
tmp <- as.matrix(wage3[,c("dlwage","dunion","dexper")])

D <- tmp[,-1]  ##D_{it}に該当する部分を取り出す
DY <- as.matrix(apply(sweep(D,1,tmp[,1],"*"),2,sum))
dhat <- solve(t(D)%*%D)%*%DY
dhat

##lm関数による推定
round(lm(dlwage~dunion+dexper,data=wage3)$coef,5)

A matrix: 2 × 1 of type dbl

dunion	0.08559373
dexper	0.06354095

(Intercept)

0

dunion

0.08559

dexper

0.06354

パッケージplmを使っても計算できる．結果は一致しているはずである．結果の見方は基本的にlm関数と同様．Fixed-effects estimatorを計算したければ，withinを指定する．indexには個体IDと時点を指定する．union, experどちらも有意だが，適合度は悪い．

result <- plm(lwage~union+exper,data=wagepan,index=c("nr","year"),model="within")
summary(result)

Oneway (individual) effect Within Model

Call:
plm(formula = lwage ~ union + exper, data = wagepan, model = "within", 
    index = c("nr", "year"))

Balanced Panel: n = 545, T = 8, N = 4360

Residuals:
     Min.   1st Qu.    Median   3rd Qu.      Max. 
-4.148119 -0.126476  0.013471  0.163406  1.505664 

Coefficients:
       Estimate Std. Error t-value  Pr(>|t|)    
union 0.0855937  0.0194323  4.4047 1.088e-05 ***
exper 0.0635409  0.0023403 27.1508 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Total Sum of Squares:    572.05
Residual Sum of Squares: 477.82
R-Squared:      0.16472
Adj. R-Squared: 0.045115
F-statistic: 375.973 on 2 and 3813 DF, p-value: < 2.22e-16

Exercise #7 (1)

1. モデル ${\tt lwage}_{it} = \delta_{1}{\tt union}_{it} + \delta_{2}{\tt exper}_{it} + u_{i} + \varepsilon_{it}$において， Pooled OLSをパッケージを利用せずに計算せよ．
2. wagepanには，1981年-1987年を表すダミー変数 d81 ~ d87も含まれている．これらを用いて，新たなモデル $$ {\tt lwage}_{it} = \delta_{1}{\tt union}_{it} + \delta_{2}{\tt exper}_{it} + \delta_{3}{\tt d82}_{it} + \delta_{4}{\tt d83}_{it}+ \cdots +\delta_{8}{\tt d87}_{it} + u_{i} + \varepsilon_{it} $$ を考え，$\delta = (\delta_{1},\dots, \delta_{9})^{T}$をfixed-effects (within) estimatorで推定せよ．
3. ダミー変数d82 ~ d87を新たに追加することのメリットを考察せよ．
4. ダミー変数d81をモデルに追加しなかった理由を考察せよ．

---

7.2 Difference-in-Differences

Settings

• $T=1$ : after treatment, $T=0$ : before treatment
• $D_{i}=1$ : treated group, $D_{i}=0$ : control, $D_{i}(t) = D_{i}I(T=t)$
• $\require{color}\textcolor{blue}{Y_{i}^{d}(t)}$ : potential outcome for $T=t$ and $D_{i}(t)=d$
• $Y_{i}(t)$ : observed outcome when $T=t$

SUTVA

• $Y_{i}(t) = D_{i}(t)Y_{i}^1(t) + (1-D_{i}(t))Y_{i}^0(t)$
  • 処置前($T=0$)においては$Y_{i}(0) = Y_{i}^0(0)$である．なぜなら，$D_{i}(0)=0$

Goal : ATT，すなわち，$E\{Y_{i}^1(1) - Y_{i}^0(1)\,|\,D_{i}=1\}$の推定

• 以後，独立なユニットを仮定 $\rightarrow$ 添字$i$は省略

Difference-in-Differences (DiD)

$$ {\rm DiD} = E\{Y(1)\,|\,D=1\}-E\{Y(0)\,|\,D=1\}-\big[E\{Y(1)\,|\,D=0\}-E\{Y(0)\,|\,D=0\}\big] $$

• 処置群におけるafeter-beforeから対照群における同一のものを引いている

DiDは平行トレンド仮定

$$ E\{Y^0(1)\,|\,D=1\} - E\{Y^0(0)\,|\,D=1\} = E\{Y^0(1)\,|\,D=0\}-E\{Y^0(0)\,|\,D=0\} $$

の下で，ATTと一致する．すなわち，$\require{color}\textcolor{blue}{{\rm DiD} = {\rm ATT}}$．よってDiDを推定すればok.

具体的にデータで計算

まずはデータを読み込む．Card and Krueger (1994)による最低賃金と雇用者数に関するパネルデータを利用する．データはhttps://www.ssc.wisc.edu/~bhansen/econometrics/ より取得し，csv形式に変換している． 解析の目的は，"最低賃金を増やすと雇用は減るか"を検証するものであり， New Jersey(NJ)は1992年に最低賃金を引き上げ，Pennsylvania(PA)はそのままである． つまり，処置群がNJ，対照群がPAである．

下記の変数のみ利用する．

• store : 店舗ID
• state : New Jerseyなら1, Pennsylvaniaなら0
• time : 最低賃金の引き上げ前なら0, 引き上げ後なら1
• empft : フルタイム雇用者数 (outcome)

## 作業ディレクトリに各自データを置いてください
ck94 <- read.csv("CK1994.csv")
ck94 <- ck94[,c("store","state","time","empft")]
head(ck94)

A data.frame: 6 × 4

	store	state	time	empft
	<int>	<int>	<int>	<chr>
1	46	0	0	30
2	49	0	0	7
3	506	0	0	3
4	56	0	0	20
5	61	0	0	6
6	62	0	0	0

変数empftには欠測値(.で表されている)があるので，その店舗を削除する．その際，同一の店舗IDが2つあることに注意(beforeとafter)．つまり，どちらかが欠測している場合は，両方とも削除しなければならない．

cck94 <- ck94[-which(ck94[,"empft"]=="."),]
ids <- which(table(cck94$store)==1) #各店舗IDが2つあるかをチェック．1つの場合の店舗IDを取り出す．
ids <- as.numeric(names(ids)) #数値にしておく
cck94 <- cck94[-which(cck94$store %in% ids),] #idsに該当する店舗の行をまとめて削除

#演算が出来るようにempftを数値型へと変換
cck94 <- transform(cck94,empft=as.numeric(empft))

7.2.1 直接的なDiDの推定

上記で定義されたDiDのnaiveな推定量は以下となる :

$$ \widehat{\rm DiD} = \frac{1}{N_{1}}\sum_{D_{i}=1}Y_{i}(1) - \frac{1}{N_{1}}\sum_{D_{i}=1}Y_{i}(0) -\bigg[ \frac{1}{N_{0}}\sum_{D_{i}=0}Y_{i}(1) - \frac{1}{N_{0}}\sum_{D_{i}=0}Y_{i}(0) \bigg] $$

これを計算するためには，StateとTimeの各組み合わせにおいてOutcomeの平均値を計算すればよい．

( tab <- cck94 %>% group_by(state,time) %>% summarize(g.empft=mean(empft)) )

`summarise()` has grouped output by 'state'. You can override using the `.groups` argument.

A grouped_df: 4 × 3

state	time	g.empft
<int>	<int>	<dbl>
0	0	10.355263
0	1	7.605263
1	0	7.743671
1	1	8.430380

(did1 <- tab$g.empft[4] - tab$g.empft[3] - (tab$g.empft[2] - tab$g.empft[1]))

3.4367088607595

7.2.2 線形回帰分析による推定

線形回帰モデルでもDiDの推定が可能．

• ${\tt State}_{i}$ : 店舗iのある地域を表すダミー変数
• ${\tt Time}_{t}$ : 時点のダミー変数． ${\tt Time}_{1}=1, {\tt Time}_{0}=0$

Potential outcomesに対する線形モデル

$$ \begin{align*} \begin{cases} Y_{i}^0(t) = \beta_{0} + \beta_{1}{\tt Comp}_{i} + \beta_{2}{\tt Time}_{t} + \varepsilon_{it} \\ Y_{i}^1(t) = Y_{i}^{0}(t) + \delta D_{i}(t) \end{cases} \end{align*} $$

ただし，$E(\varepsilon_{it}\,|\,D_{i}(1))=0$とする．このとき，SUTVAから，

$$ Y_{i}(t) = \beta_{0} + \beta_{1}{\tt State}_{i} + \beta_{2}{\tt Time}_{t} + \delta D_{i}(t) + \varepsilon_{it} \tag{1} $$

であり，$\require{color}\textcolor{red}{{\rm DiD} = \delta}$である(実際に計算してみよ)． よってlm関数で簡単にDiDの推定ができる．

##state*time=Di(t)であることに注意
summary(lm(empft~state+time+state*time,data=cck94))

Call:
lm(formula = empft ~ state + time + state * time, data = cck94)

Residuals:
   Min     1Q Median     3Q    Max 
-10.36  -6.43  -2.43   3.57  52.26 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  10.3553     0.9545  10.849   <2e-16 ***
state        -2.6116     1.0631  -2.457   0.0142 *  
time         -2.7500     1.3498  -2.037   0.0420 *  
state:time    3.4367     1.5034   2.286   0.0225 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.321 on 780 degrees of freedom
Multiple R-squared:  0.00844,	Adjusted R-squared:  0.004626 
F-statistic: 2.213 on 3 and 780 DF,  p-value: 0.08523

Exercise #7 (2)

式(1)を適当に変形すると，下記のようなfixed-effects (within)モデルへと変形できる．

$$ Y_{it} = \delta^{T}D_{it} + u_{i} + \varepsilon_{it} $$

ここで$\delta$および$D_{it}$はベクトルになり得る．式(1)を正しく変形し，次の問いに答えよ．

1. 関数plmを用いてDiDの推定値を求めよ．
2. 関数plmを用いずにDiDの推定値を計算せよ．

---

7.3 Synthetic Control Method

パッケージtidysynthの読み込み．California's proportion 99のデータを利用する．

library(tidysynth)
data(smoking) #California's proportion 99

smokingはロング型のデータなので，ワイド型のデータに変換する．変換にはtidyrライブラリ内のpivot_wider関数を利用する．names_fromにはワイドに展開したい変数を，values_fromには展開したい値をそれぞれ指定する．

wsmoking <- smoking[,c("state","year","cigsale")] %>% pivot_wider(names_from=year,values_from=cigsale)

使いやすいようにdata.frame型に変形．

data <- data.frame(wsmoking)
colnames(data) <- colnames(wsmoking)
treat <- which(data[,"state"]=="California")  #treatment group
t0 <- which(colnames(data[,-1])=="1988")  #T0

y10 <- as.vector(data[treat,-1])[1:t0]  #Carifornia before treatment
y00 <- as.matrix(data[-treat,-1])[,1:t0] #Others before treatment
y11 <- as.vector(data[treat,-1])[(t0+1):31] #Carifornia after treatment
y01 <- as.matrix(data[-treat,-1])[,(t0+1):31] #Others after treatment

7.3.1 Ridge回帰によるSynthetic controlの計算

library(glmnet)
result <- glmnet(x=t(y00),y=t(y10),alpha=0,lambda=0.01) #lambda=0.01でridge回帰を実行
#coef(result)

全ての期間において，Synthetic control $\sum_{i=2}^{N}w_{i}^{*}Y_{i}(t)$を計算．

synth <- predict(result,newx=t(as.matrix(data[-treat,-1])))

全期間におけるSynthetic controlと実際のCaliforniaの売上データを比較する．

origin.ca <- as.vector(data[treat,-1]); origin.ca <- t(origin.ca)

plot(1970:2000,origin.ca,ylim=c(38,130),xlab="",ylab="tabaco sales",type="l",lty=1,lwd=2)
par(new=T)
plot(1970:2000,synth,ylim=c(38,130),xlab="",ylab="tabaco sales",type="l",lty=1,col="blue",lwd=2)
abline(v=1988)

7.3.2 凸包制約付き最適化によるSynthetic controlの計算

$$ \min\sum_{k=1}^{K}v_{k}\bigg\{X_{1k} - \sum_{i=2}^{N}w_{i}X_{ik}\bigg\}^{2} \quad \mbox{subject to}\quad w_{2},\dots, w_{N} \ge 0\quad\mbox{and}\quad \sum_{i=2}^{N}w_{i}=1 $$

簡単のため$v_{k}=1$とする．solve.QP関数で制約付き最適化問題を解く．

library(quadprog)

Dmat <- y00 %*% t(y00) + 0.01*diag(38)  #Positive definiteになるように補正
dvec <- y00 %*% t(y10)
Amat <- t(rbind(rep(1,38),diag(38)))
bvec <- as.vector(c(1,rep(0,38)))

result <- solve.QP(Dmat,dvec,Amat,bvec,meq=1)
what <- data.frame(round(result$solution,4)) #解
colnames(what) <- "weight"; rownames(what) <- data$state[-treat] #名前の変更
what

A data.frame: 38 × 1

	weight
	<dbl>
Rhode Island	0.0000
Tennessee	0.0000
Indiana	0.0000
Nevada	0.2049
Louisiana	0.0000
Oklahoma	0.0000
New Hampshire	0.0454
North Dakota	0.0000
Arkansas	0.0000
Virginia	0.0000
Illinois	0.0000
South Dakota	0.0000
Utah	0.3939
Georgia	0.0000
Mississippi	0.0000
Colorado	0.0148
Minnesota	0.0000
Texas	0.0000
Kentucky	0.0000
Maine	0.0000
North Carolina	0.0000
Montana	0.2318
Vermont	0.0000
Iowa	0.0000
Connecticut	0.1091
Kansas	0.0000
Delaware	0.0000
Wisconsin	0.0000
Idaho	0.0000
New Mexico	0.0000
West Virginia	0.0000
Pennsylvania	0.0000
South Carolina	0.0000
Ohio	0.0000
Nebraska	0.0000
Missouri	0.0000
Alabama	0.0000
Wyoming	0.0000

Synthetic controlの作成には，Nevada, New Hampshire, Utah, Colorado, Montana, Connecticutのみが効いていることがわかる．おそらくこれらの州が，処置前期間においてCaliforniaのたばこ売上と近い．

同様にSynthetic controlをプロットしてみる.今度は処置前期間において過剰適合していないし，処置後においてはっきりと差が見られる．

synth2 <- t(as.matrix(data[-treat,-1])) %*% what$weight

plot(1970:2000,origin.ca,ylim=c(38,130),xlab="",ylab="tabaco sales",type="l",lty=1,lwd=2)
par(new=T)
plot(1970:2000,synth2,ylim=c(38,130),xlab="",ylab="tabaco sales",type="l",lty=1,col="red",lwd=2)
abline(v=1988)

7.3.3 Fisherの正確検定によるp値の計算

上記の処置後期間における差が偶然でないかどうかを，Fisherの正確検定を用いて検証する．

p値は次のようなステップで計算される．

1. $N=39$州のそれぞれを仮に処置群だと想定する．
2. 各データセット($N=39$)についてSynthetic control，すなわち，$w_{i}$の推定を行う.
3. 処置前期間におけるMSPEをそれぞれ($N=39$)について計算． $$ MSPE_{\tt \,before} = \frac{1}{T_{0}}\sum_{t=1}^{T_{0}} \big\{Y_{1}(t) - \sum_{i=2}^{N}\hat{w}_{i}Y_{i}(t)\big\}^{2} $$
4. 処置後期間においても同様にMSPEを計算する． $$ MSPE_{\tt \,after} = \frac{1}{T-T_{0}}\sum_{t=T_{0}+1}^{T} \big\{Y_{1}(t) - \sum_{i=2}^{N}\hat{w}_{i}Y_{i}(t)\big\}^{2} $$
5. 比$MSPE_{\tt \,after}/MSPE_{\tt \,before}$を計算する．
6. 計算した比を降順にソートし，CAの順位$r$を求める．
7. 最後にp値 : $p=\frac{r}{N}$を計算する．

Exercise #7 (3)

上記をCalifornia's proportion 99のデータで実行し，CAにおけるp値を算出せよ．