Instrumental Variable Regression(도구변수를 활용한 회귀분석)3

이전 포스팅(direction-f.tistory.com/57, direction-f.tistory.com/58)에서 도구변수의 개념 및 일반화된 도구변수를 활용한 회귀분석과 도구변수의 Validity를 검증할 수 있는 방안에 대해서 간략히 정리하였습니다.

이번 포스팅에서는, 예제를 통해 도구변수를 활용한 회귀분석을 수행하고 타당성 검증을 수행해보겠습니다.

우리는 위와 같은 모형을 추정해보겠습니다.

먼저 필요한 데이터를 만들고 모형을 추정해보겠습니다. 여기서 내생변수는 pricediff(코드상 endog)입니다.

import pandas as pd
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
import seaborn as sns
from linearmodels.iv import IV2SLS

## Application
c1985 = cigaret.loc[cigaret["year"]==1985]
c1995 = cigaret.loc[cigaret["year"]==1995]

packdiff = np.log(c1995["packs"].values)-np.log(c1985["packs"].values)
pricediff = np.log(c1995["rprice"].values)-np.log(c1985["rprice"].values)
incomediff = np.log(c1995["rincome"].values)-np.log(c1985["rincome"].values)

salestaxdiff = c1995["salestax"].values-c1985["salestax"].values
cigataxdiff = c1995["cigtax"].values-c1985["cigtax"].values

## Instrumental variable
dependent = packdiff
exog = sm.add_constant(incomediff)
endo = pricediff

#iv-1
instrument1 = salestaxdiff
cig_ivreg_diff1 = IV2SLS(dependent, exog, endo, instrument1).fit()
cig_ivreg_diff1 

'''
==============================================================================
            Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI
------------------------------------------------------------------------------
exog.0        -0.1180     0.0661    -1.7859     0.0741     -0.2474      0.0115
exog.1         0.5260     0.3287     1.6001     0.1096     -0.1183      1.1702
endog         -0.9380     0.2009    -4.6688     0.0000     -1.3318     -0.5442
==============================================================================
'''

#iv-2
instrument2 = cigataxdiff
cig_ivreg_diff2 = IV2SLS(dependent, exog, endo, instrument2).fit()
cig_ivreg_diff2 

'''
==============================================================================
            Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI
------------------------------------------------------------------------------
exog.0        -0.0170     0.0651    -0.2620     0.7934     -0.1446      0.1105
exog.1         0.4281     0.2892     1.4803     0.1388     -0.1387      0.9950
endog         -1.3425     0.2214    -6.0638     0.0000     -1.7764     -0.9086
==============================================================================
'''

#iv-3
instrument3 = np.column_stack((salestaxdiff,cigataxdiff))
cig_ivreg_diff3 = IV2SLS(dependent, exog, endo, instrument3).fit()
cig_ivreg_diff3 

'''
==============================================================================
            Parameter  Std. Err.     T-stat    P-value    Lower CI    Upper CI
------------------------------------------------------------------------------
exog.0        -0.0520     0.0605    -0.8595     0.3901     -0.1706      0.0666
exog.1         0.4620     0.2995     1.5426     0.1229     -0.1250      1.0491
endog         -1.2024     0.1907    -6.3056     0.0000     -1.5761     -0.8287
==============================================================================
'''

도구변수의 타당성을 검증하기 위해 내생변수(pricediff)와 도구변수간 회귀분석을 수행하겠습니다. 이 검증결과 하나 이상의 도구변수가 내생변수를 설명하는데 관련성을 가짐을 얻을 수 있어야 합니다. 어떠한 도구변수도 내생변수를 설명하는데 유의성이 없다면 잘못된 도구변수를 활용했다는 것을 알 수 있습니다.

# first-stage regressions
temp_df = pd.DataFrame(np.column_stack((pricediff, salestaxdiff, incomediff, cigataxdiff)), columns =["pricediff", "salestaxdiff", "incomediff", "cigataxdiff"])
mod_relevance1 = smf.ols("pricediff~salestaxdiff+incomediff", data=temp_df).fit()
mod_relevance2 = smf.ols("pricediff~cigataxdiff+incomediff", data=temp_df).fit()
mod_relevance3 = smf.ols("pricediff~cigataxdiff+incomediff+salestaxdiff", data=temp_df).fit()

## F-test
from statsmodels.stats.anova import anova_lm
base = smf.ols("pricediff~incomediff", data=temp_df).fit()

# check instrument relevance for model (1)
res=anova_lm(base,mod_relevance1)
print(res)
'''
   df_resid       ssr  df_diff   ss_diff          F        Pr(>F)
0      46.0  0.366762      0.0       NaN        NaN           NaN
1      45.0  0.180550      1.0  0.186212  46.411287  1.933476e-08
'''

# check instrument relevance for model (1)
res2=anova_lm(base,mod_relevance2)
print(res2)
'''
   df_resid       ssr  df_diff   ss_diff          F        Pr(>F)
0      46.0  0.366762      0.0       NaN        NaN           NaN
1      45.0  0.119190      1.0  0.247573  93.470784  1.481292e-12
'''

# check instrument relevance for model (3)
res3=anova_lm(base,mod_relevance3)
print(res3)
'''
   df_resid       ssr  df_diff   ss_diff          F        Pr(>F)
0      46.0  0.366762      0.0       NaN        NaN           NaN
1      44.0  0.082627      2.0  0.284135  75.652583  5.757803e-15
'''

위의 결과를 보면 최소한 하나 이상의 도구변수가 내생변수를 설명하는데 유의함을 알 수 있습니다. 

추가적으로, 도구변수를 활용하여 추정한 모형의 Residual과 실제 도구변수가 관계가 없음을 검증해보도록 하겠습니다.(관계가 없어야 합니다.)

앞선 포스팅에서 정리했듯이 이 때는 F-통계량 대신 J-통계량을 활용하겠습니다. 이는 $\chi2$ 통계량이며, $J=mF$입니다. 이때 $m$은 도구변수의 수입니다.

temp_df["resid"]=cig_ivreg_diff3.resids
cig_iv_or = smf.ols("resid~incomediff+salestaxdiff+cigataxdiff",data=temp_df).fit()
base= smf.ols("resid~incomediff", data=temp_df).fit()

## F-statistic
res=anova_lm(base, cig_iv_or)

## J-statistic
chi_squared_statistic=res.F[1]*2

## chi-square test
from scipy.stats import chi2
dof = 2-1
chi_ = chi2(dof)

p_val=1-chi_.cdf(chi_squared_statistic)
print(p_val) ##0.02636406040208783

J-test 결과를 보면, p-value가 0.05 이하로 도구변수와 Residual와 관계성을 가짐을 알 수 있습니다. 이는 salsestaxidff 또는 cigataxdiff가 유의한 도구변수가 아님을 뜻할 수 있습니다.