In [ ]:
# 必要なライブラリをインストール
!pip install japanize_matplotlib
In [ ]:
# ライブラリを起動
from dateutil.relativedelta import relativedelta
import japanize_matplotlib
import math as math
import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
import random
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLS
from statsmodels.tsa.api import VAR
from sklearn.linear_model import LinearRegression
from regex import R
from datetime import datetime as dt
import datetime
import os
from matplotlib import pyplot as plt
import matplotlib.dates as mdates
import matplotlib.ticker as ticker
from matplotlib.pylab import rcParams
rcParams["figure.figsize"] = 15, 6
# Colabでファイルを読み込むために、Google Driveをマウント
from google.colab import drive
drive.mount('/content/drive')
Mounted at /content/drive
In [ ]:
# データの読み込み
DATASET_VAR = pd.read_excel('/content/drive/My Drive/Data/data_VAR8.xlsx',
sheet_name='data_VAR', header=[0], index_col=0)
DATASET_IV = pd.read_excel('/content/drive/My Drive/Data/data_VAR8.xlsx',
sheet_name='data_IV', header=[0], index_col=0)
class structtype():
def __init__(self):
self.data = []
def psvar(VARresids, proxy, policyVarPos, nLag):
if np.sum(np.isnan(proxy)) > 0 and len(proxy) == len(VARresids):
index_notnan = ~np.isnan(proxy[:, 0])
proxy = proxy[index_notnan, :]
VARresids = VARresids[index_notnan, :]
VARresids = np.matrix(VARresids)
proxy = np.matrix(proxy)
T, N = VARresids.shape
t, n = proxy.shape
rows_iP = list(range(N))
rows_iP.remove(policyVarPos)
rows_iP = [policyVarPos] + rows_iP
rfResids = VARresids[:, rows_iP]
IV1X = np.concatenate([np.ones([t, 1]), proxy], axis=1)
IV1_coef = np.linalg.lstsq(IV1X, rfResids[:, 0], rcond=None)[0]
IV1_resids = rfResids[:, 0] - IV1X @ IV1_coef
IV_uphat = IV1X @ IV1_coef
IV2X = np.concatenate([np.ones([t, 1]), IV_uphat], axis=1)
IV2_coef = np.linalg.lstsq(IV2X, rfResids[:, 1:], rcond=None)[0]
IV2_resids = rfResids[:, 1:] - IV2X @ IV2_coef
s21is11 = IV2_coef[1:, :].T
SigmaU = np.cov(rfResids.T)
sigma11 = SigmaU[0:1, 0:1]
sigma12 = SigmaU[0:1, 1:]
sigma21 = SigmaU[1:, 0:1]
sigma22 = SigmaU[1:, 1:]
Qmat = np.dot(s21is11*sigma11, s21is11.T) - (
np.dot(sigma21, s21is11.T) + np.dot(s21is11, sigma21.T)) + sigma22
Qinv = np.linalg.solve(Qmat, (sigma21 - s21is11 * sigma11))
s12s12p = np.dot((sigma21- s21is11*sigma11).T, Qinv)
s11s11p = sigma11 - s12s12p
s11 = math.sqrt(s11s11p)
s1 = np.empty([N, 1])
s1[0, 0] = s11
s1[1:, 0:1] = s21is11*s11
IV1_dev = rfResids[:, 0] - np.mean(rfResids[:, 0])
SST_m = IV1_dev.T @ IV1_dev
SSE_m = IV1_resids[:, 0].T @ IV1_resids[:, 0]
Fstat = ((SST_m - SSE_m) / n) / (SSE_m / (t - n - 1))
R2 = 1 - SSE_m / SST_m
R2adj = R2 - (n / (t - n - 1)) * (1 - R2)
SS_m = np.zeros([n + 1, n + 1])
for ii in range(t):
SS_m = SS_m + (1 / t * IV1X[ii, :].T @ IV1X[ii, :] * (
IV1_resids[ii, 0] ** 2))
Avarb_m = np.linalg.inv(1 / t * IV1X.T @ IV1X) * SS_m * np.linalg.inv(
1 / t * IV1X.T @ IV1X)
IV1_robSE = np.sqrt(Avarb_m[1, 1])
IV1_robtstat = IV1_coef[1] / IV1_robSE
RR_m = np.concatenate([np.zeros([n, 1]), np.eye(n)], axis=1)
WW_m = t * (RR_m@IV1_coef[:, 0]).T @ np.linalg.inv(
RR_m @ Avarb_m @ RR_m.T) @ (RR_m @ IV1_coef[:, 0])
Fstat_rob = WW_m / n
res = structtype
res.B = s1[1:, :]
res.B = np.insert(res.B, policyVarPos, s1[0, 0], axis=0)
res.Coef = np.asarray(IV1_coef).reshape(-1).tolist()
res.Coef_const = np.asarray(IV1_coef[0, 0]).reshape(-1).tolist()
res.Obs = t
res.R2 = R2[0, 0]
res.R2adj = R2adj[0, 0]
res.Fstat = Fstat[0, 0]
res.Fstat_rob = Fstat_rob[0, 0]
Bzero = np.zeros([N, N - 1])
Bzero = np.insert(Bzero, policyVarPos, res.B.T, axis=1)
model_ols = OLS(rfResids[:,0], IV1X)
result = model_ols.fit(cov_type = 'HC0')
res.SE = result.bse
res.tstat = result.tvalues.tolist()
res.pvalue = result.pvalues.tolist()
result.fvalue
result.rsquared
result.rsquared_adj
return res, Bzero
VAR_varname = ['IP', 'cCPItaxSA', 'JGB1Y', 'TOPIXMA', 'CBSm']
VAR_figname = ['鉱工業生産', '消費者物価指数', '1年物国債利回り', '東証株価指数',
'社債スプレッド']
IV_varname = ['TARGET']
MPvar = 'JGB1Y'
VAR_smpl_from = '1990-01-01'
VAR_smpl_until = '2020-01-01'
IV_smpl_from = '1992-07-01'
IV_smpl_until = '2020-01-01'
iScheme = 'IV'
nLag = 12
nHor = 48
bandsCoverage = 0.9
shockSize = 1
def Gmatrices(AL, C, nLag, nHor, n):
J = np.hstack([I, np.zeros([n, (nLag-1)*n])])
Alut = np.vstack(
[AL, np.hstack([np.eye(n*(nLag-1)), np.zeros([n*(nLag-1), n])])])
AJ = np.zeros([n*nLag, n, nHor+1])
for j in range(nHor+1):
if j == 0:
AJ[:, :, j] = np.eye(Alut.shape[0]) @ J.T
else:
AJ[:, :, j] = np.linalg.matrix_power(Alut, j) @ J.T
JAp = np.reshape(AJ, [n*nLag, n*(nHor+1)], order='F').T
AJaux = np.zeros([JAp.shape[0]*n, JAp.shape[1]*n, nHor+1])
Caux = np.reshape(
np.hstack([I, C[:, :nHor*n]]), [n, n, nHor+1], order='F')
for i in range(nHor+1):
AJaux[((n**2)*i):, :, i] = np.kron(JAp[:n*(nHor+1-i), :],
Caux[:, :, i])
Gaux = np.transpose(np.reshape(
np.sum(AJaux, axis=2).T, [(n**2)*nLag, n**2, nHor+1], order='F'),
[1, 0, 2])
G = np.zeros([Gaux.shape[0], Gaux.shape[1], Gaux.shape[2]+1])
G[:, :, 1:] = Gaux
return Gaux
try:
VAR_smpl_from_index = DATASET_VAR.index.get_loc(VAR_smpl_from)
VAR_smpl_until_index = DATASET_VAR.index.get_loc(VAR_smpl_until)
except KeyError:
print('ERROR: CHECK THE VAR SAMPLE PERIOD YOU SET.')
if np.count_nonzero(np.isin(DATASET_VAR.columns, VAR_varname)) != len(
VAR_varname):
print('ERROR: CHECK THE VAR VARIABLE NAMES YOU SET.')
Ydf = DATASET_VAR.loc[VAR_smpl_from:VAR_smpl_until, VAR_varname]
print('VAR VARIABLE:', *Ydf.columns, sep=' ')
print('VAR SAMPLE: ', dt.date(Ydf.index[0]), '-', dt.date(Ydf.index[-1]))
if iScheme == 'IV':
IVstart = DATASET_IV[IV_varname].dropna().index[0]
if dt.strptime(IV_smpl_from, '%Y-%m-%d') < IVstart:
IV_smpl_from = dt.strftime(IVstart, '%Y-%m-%d')
try:
VAR_smpl_from_IV_index = DATASET_VAR.index.get_loc(IV_smpl_from)
VAR_smpl_until_IV_index = DATASET_VAR.index.get_loc(IV_smpl_until)
except KeyError:
print('ERROR: The IV sample is longer than the VAR sample period.')
if VAR_smpl_from_IV_index < VAR_smpl_from_index + nLag:
VAR_smpl_from_IV_index = VAR_smpl_from_index + nLag
IV_smpl_from = DATASET_VAR.index[VAR_smpl_from_IV_index]
if VAR_smpl_until_IV_index > VAR_smpl_until_index:
IV_smpl_until = VAR_smpl_until
VAR_smpl_until_IV_index = VAR_smpl_until_index
IVindex = np.empty([2, 1], dtype='int')
IVindex[0] = VAR_smpl_from_IV_index - (VAR_smpl_from_index + nLag)
IVindex[1] = VAR_smpl_until_index - VAR_smpl_until_IV_index
try:
IV_smpl_from_index = DATASET_IV.index.get_loc(IV_smpl_from)
IV_smpl_until_index = DATASET_IV.index.get_loc(IV_smpl_until)
except KeyError:
print('ERROR: CHECK THE IV SAMPLE PERIOD YOU SET.')
if np.count_nonzero(np.isin(DATASET_IV.columns, IV_varname)) != len(
IV_varname):
print('ERROR: CHECK THE INSTRUMENT NAME YOU SET.')
IVdf = DATASET_IV.loc[IV_smpl_from:IV_smpl_until, IV_varname]
print('IV VARIABLE:', *IVdf.columns, sep=' ')
print('IV SAMPLE: ', dt.date(IVdf.index[0]), '-', dt.date(IVdf.index[-1]))
elif iScheme == 'CHOL':
print('IDENTIFICATION: Cholesky decomposition')
if Ydf.iloc[0, :].isnull().sum() > 0:
print('*** ERROR: CHECK THE BEGINNING OF THE VAR SAMPLE PERIOD ***')
if Ydf.iloc[-1, :].isnull().sum() > 0:
print('*** ERROR: CHECK THE END OF THE VAR SAMPLE PERIOD ***')
Y = Ydf.values
if iScheme == 'IV':
IV = IVdf.values
MPshockIndex = VAR_varname.index(MPvar)
T, n = Y.shape
shockS = np.zeros([n, 1])
sLevel = abs(stats.norm.ppf((1-bandsCoverage)/2))
IRF = np.empty([n, nHor+1])
IRF.fill(np.nan)
Yh = Y[nLag:, :]
nT, nY = Yh.shape
YhLag = np.empty([T-nLag, n*nLag])
YhLag.fill(np.nan)
for j in range(nLag):
YhLag[:, n*j:n*(j+1)] = Y[nLag-j-1:-j-1, :].copy()
YhprojSet = np.concatenate([np.ones([nT, 1]), YhLag], axis=1)
projCoeffs = np.linalg.lstsq(YhprojSet, Yh, rcond=None)[0]
irfCoeffs = projCoeffs[1:n+1, :]
u = Yh - YhprojSet @ projCoeffs
SigmaU = np.cov(u.T)
if iScheme == 'CHOL':
Bzero = np.linalg.cholesky(SigmaU)
elif iScheme == 'IV':
if IVindex[1, 0] == 0:
FirstStage, Bzero = psvar(
u[IVindex[0,0]:, :], IV, MPshockIndex, nLag)
else:
FirstStage, Bzero = psvar(
u[IVindex[0, 0]:-IVindex[1, 0], :], IV, MPshockIndex, nLag)
print('F-stat: ', round(FirstStage.Fstat, 3),
', robust F-stat:', round(FirstStage.Fstat_rob, 3))
print('R2:', round(FirstStage.R2, 3))
print('Adj R2:', round(FirstStage.R2adj, 3))
print('SE:', [round(i, 3) for i in FirstStage.SE])
print('Coef:', [round(i, 3) for i in FirstStage.Coef])
print('t-stat:', [round(i, 3) for i in FirstStage.tstat])
print('p-value:', [round(i, 3) for i in FirstStage.pvalue])
FirstStage_IV_Fstat = FirstStage.Fstat.copy()
FirstStage_IV_Fstat_rob = FirstStage.Fstat_rob.copy()
FirstStage_IV_R2 = FirstStage.R2.copy()
FirstStage_IV_R2adj = FirstStage.R2adj.copy()
FirstStage_IV_SE = FirstStage.SE.copy()
FirstStage_IV_Coef = FirstStage.Coef.copy()
FirstStage_IV_tstat = FirstStage.tstat.copy()
FirstStage_IV_pvalue = FirstStage.pvalue.copy()
if shockSize == False:
shockS[MPshockIndex] = 1
else:
shockS[MPshockIndex] = shockSize / Bzero[MPshockIndex, MPshockIndex]
IRF[:, 0:1] = Bzero @ shockS
A = np.zeros([n*nLag, n*nLag])
A[0:n, :] = projCoeffs[1:, :].T
A[n:, 0:n*(nLag-1)] = np.eye(n*(nLag-1))
for h in range(nHor):
Ah = np.linalg.matrix_power(A, h+1)
IRF[:, h+1:h+2] = Ah[0:n, 0:n] @ Bzero @ shockS
m = YhprojSet.shape[1] - (n*nLag)
XSVARp = YhprojSet
if iScheme == 'IV':
k = IV.shape[1] + 1
IVm = np.empty([Y.shape[0], k])
IVm.fill(np.nan)
if IVindex[1, 0] == 0:
IVm[IVindex[0, 0]+nLag:, 0] = 1
IVm[IVindex[0, 0]+nLag:, 1:] = IV
else:
IVm[IVindex[0, 0]+nLag:-IVindex[1, 0], 0] = 1
IVm[IVindex[0, 0]+nLag:-IVindex[1, 0], 1:] = IV
IVm = IVm[nLag:, :]
IVm = np.nan_to_num(IVm)
else:
k=1
IVm = np.zeros([349, 1])
matagg = np.hstack([XSVARp, u, IVm]).T
T1aux = u.shape[0]
T2aux = matagg.shape[0]
etaaux = np.reshape(u.T, [nY, 1, T1aux], order='F')
mataggaux = np.transpose(np.reshape(matagg, [T2aux, 1, T1aux]), [1, 0, 2])
auxeta = np.multiply(etaaux, mataggaux) - \
np.mean(np.multiply(etaaux, mataggaux), 2)[:, :, np.newaxis]
vecAss1 = np.reshape(
auxeta, [(nY*m)+(nLag*(nY**2))+nY**2+(nY*k), 1, T1aux], order='F')
AuxHAC1 = vecAss1
AuxHAC2 = np.reshape(
AuxHAC1, [vecAss1.shape[0], vecAss1.shape[2]], order='F').T
NWlags = int(np.floor(4*((nT/100) ** (2/9))))
Sigma0 = (1/(AuxHAC2.shape[0]))*(AuxHAC2.T@AuxHAC2)
SigmaNW = Sigma0.copy()
for i in range(NWlags):
Sigma_cov_i = (1/AuxHAC2.shape[0]) * (AuxHAC2[:-i-1, :].T@AuxHAC2[1+i:, :])
SigmaNW = SigmaNW + (1-(i+1)/(NWlags+1)) * (Sigma_cov_i+Sigma_cov_i.T)
AuxHAC3 = SigmaNW
WhatAss1 = AuxHAC3
I = np.eye(nY)
V = np.kron(I[0, :], I)
for i_vars in range(1, n):
V = np.vstack([V, np.kron(I[i_vars, :], I[i_vars:, :])])
Q1 = XSVARp.T @ XSVARp / T1aux
Q2 = IVm.T @ XSVARp / T1aux
invQ1 = np.linalg.inv(Q1)
Shat_a1 = np.kron(np.hstack([np.zeros([n*nLag, m]), np.eye(n*nLag)])@invQ1, I)
Shat_a2 = np.zeros([(n**2)*nLag, n**2+(k*n)])
Shat_b = np.hstack([np.zeros([int(n*(n+1)/2), ((n**2)*nLag)+n *
m]), V, np.zeros([int(n*(n+1)/2), k*n])])
Shat_c = np.hstack(
[- np.kron(Q2@invQ1, I), np.zeros([k*n, n**2]), np.eye(k*n)])
Shat = np.vstack([np.hstack([Shat_a1, Shat_a2]), Shat_b, Shat_c])
WHataux = Shat @ WhatAss1 @ Shat.T
WHat11 = WHataux[:(n**2)*nLag, :(n**2)*nLag]
WHat12 = WHataux[:(n**2)*nLag, ((n**2)*nLag)+int(n*(n+1)/2):]
WHat22 = WHataux[((n**2)*nLag)+int(n*(n+1)/2):,
((n**2)*nLag)+int(n*(n+1)/2):]
WHatSigma = WHataux[(n**2)*nLag:((n**2)*nLag)+int(n*(n+1)/2),
(n**2)*nLag:((n**2)*nLag)+int(n*(n+1)/2)]
What = np.vstack([np.hstack([WHat11, WHat12]),
np.hstack([WHat12.T, WHat22])])
AL = projCoeffs[1:, :].T
vecAL = np.reshape(AL, [n, n, nLag], order='F')
vecALrevT = np.zeros([n, n, nHor])
for i in range(nHor):
if i < (nHor-nLag):
vecALrevT[:, :, i] = np.zeros([n, n])
else:
vecALrevT[:, :, i] = vecAL[:, :, nHor-i-1].T
vecALrevT = np.reshape(vecALrevT, [n, n*(nHor)], order='F')
MARepC = np.tile(vecAL[:, :, 0], [1, nHor])
for i in range(nHor-1):
MARepC[:, (n*(i+1)):(n*(i+2))] = np.hstack([I, MARepC[:, :n*(i+1)]]
) @ vecALrevT[:, (
nHor*n-(n*(i+2))):].T
Caux = np.hstack([I, MARepC])
C = np.reshape(Caux, [n, n, nHor+1], order='F')
G = Gmatrices(AL, Caux, nLag, nHor, n)
scale = 1
clevel = bandsCoverage
critval = sp.stats.norm.ppf(1-(1-clevel)/2, 0, 1) ** 2
lambdahat = np.zeros([n, nHor+1])
DmethodVar = np.zeros([n, nHor+1])
Dmethodlbound = np.zeros([n, nHor+1])
Dmethodubound = np.zeros([n, nHor+1])
if iScheme == 'IV':
Gamma_m = np.zeros([n, k])
for i in range(k):
Gamma_m[:,i] = u.T @ IVm[:,i] / nT
BzeroMP = Bzero[:, MPshockIndex][:, np.newaxis]
if shockSize == False:
Bzero00 = 1
else:
Bzero00 = shockSize*BzeroMP[MPshockIndex, 0]
for ih in range(nHor+1):
for ivar in range(n):
lambdahat[ivar, ih] = scale * \
I[:, ivar].T @ C[:, :, ih] @ BzeroMP / Bzero00
d1 = scale * np.kron(BzeroMP.T, I[:, ivar].T) @ G[:, :, ih]
sm = np.kron(Gamma_m[ivar,:], (I[:, ivar].T @ C[
:, :, ih] - lambdahat[ivar, ih] * I[:, MPshockIndex].T))
d = np.vstack([d1.T, sm[:, np.newaxis]])
DmethodVar[ivar, ih] = d.T @ What @ d
Dmethodlbound[ivar, ih] = lambdahat[ivar, ih] - \
((critval/nT)**.5)*(DmethodVar[ivar, ih]**.5)/abs(Bzero00)
Dmethodubound[ivar, ih] = lambdahat[ivar, ih] + \
((critval/nT)**.5)*(DmethodVar[ivar, ih]**.5)/abs(Bzero00)
else:
Bzero00 = Bzero[MPshockIndex, MPshockIndex]
w = np.arange(nY*nY).reshape((nY, nY), order='F').T.ravel(order='F')
Kkk = np.eye(nY*nY)[w, :]
H = V.T @ np.linalg.inv(V @ (np.eye(nY**2)+Kkk) @ np.kron(Bzero,I) @ V.T)
for ih in range(nHor+1):
lambdahat[:, ih] = scale * I @ C[:, :, ih] @ Bzero[
:, MPshockIndex] / Bzero00
lambdahat_all_ih = scale * I @ C[:, :, ih] @ Bzero
if ih == 0:
d1 = np.zeros([nY**2,nLag*nY**2])
else:
d1 = scale * np.kron(Bzero.T, I) @ G[:, :, ih]
d2 = scale * np.kron(I, lambdahat_all_ih) @ H
DmethodVar_h = d1 @ WHat11 @ d1.T + d2 @ WHatSigma @ d2.T
DmethodVar[:, ih] = np.diag(DmethodVar_h)[
(MPshockIndex-1)*nY:MPshockIndex*nY]
Dmethodlbound[:, ih] = lambdahat[:, ih] - \
((critval/nT)**.5)*(DmethodVar[:, ih]**.5) / Bzero00
Dmethodubound[:, ih] = lambdahat[:, ih] + \
((critval/nT)**.5)*(DmethodVar[:, ih]**.5) / Bzero00
IRFl = Dmethodlbound.copy()
IRFh = Dmethodubound.copy()
IRF = lambdahat
if iScheme == 'CHOL':
figcolor = 'r'
elif iScheme == 'IV':
figcolor = 'b'
if not 'VAR_figname' in locals():
VAR_figname = VAR_varname
fig_hrzn = np.arange(0, nHor+1)
nvars = len(VAR_varname)
fig = plt.figure(figsize=(nvars*3, 8))
ylim=([-10, 6], [-1.4, 0.4], [0, 1.8], [-50, 10], [-0.4, 0.6])
for i in range(nvars):
ax = fig.add_subplot(2, 3, i+1)
ax.set_title(VAR_figname[i], fontsize=14)
ax.axhline(y=0, xmin=0, xmax=1, color='black', linewidth=0.75)
ax.plot(fig_hrzn, IRF[i, :], color='black', linewidth=2)
ax.plot(fig_hrzn, IRFl[i, :], color='black', linewidth=1,
linestyle='dashed')
ax.plot(fig_hrzn, IRFh[i, :], color='black', linewidth=1,
linestyle='dashed')
ax.set_xlim([0, nHor])
ax.set_xlabel('月', fontsize=14)
ax.set_xticks(range(0, nHor+1, 6))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_ylim(ylim[i])
fig.tight_layout()
SaveDirName = 'Figure/VAR'
if not os.path.exists('Figure'):
os.makedirs('Figure')
if not os.path.exists(SaveDirName):
os.makedirs(SaveDirName)
if iScheme == 'CHOL':
FigName = 'irsSVAR_JP_Chol'
elif iScheme == 'IV':
FigName = 'irsSVAR_JP_IV'
SaveFigName = SaveDirName + '/' + FigName
with open(SaveFigName+'.txt', 'w') as f:
f.write('LAG ORDER: ' + str(nLag) + '\n')
f.write('POLICY INDICATOR: ' + MPvar + '\n')
f.write('VAR VARIABLE: ' + str(Ydf.columns.values) + '\n')
f.write('VAR SAMPLE: ' + str(dt.date(Ydf.index[0])) + ' - ' + str(
dt.date(Ydf.index[-1])) + '\n')
if iScheme == 'IV':
f.write('IV VARIABLE: ' + str(IVdf.columns.values) + '\n')
f.write('IV SAMPLE: ' + str(dt.date(IVdf.index[0])) + ' - ' + str(
dt.date(IVdf.index[-1])) + '\n')
f.write('F-stat: ' + str(round(FirstStage_IV_Fstat, 3)
) + ', robust F-stat: ' + str(round(FirstStage_IV_Fstat_rob, 3)) + '\n')
f.write('R2: ' + str(round(FirstStage_IV_R2, 3)) + '\n')
f.write('Adj R2: ' + str(round(FirstStage_IV_R2adj, 3)) + '\n')
f.write('SE: ' + str([round(i, 3) for i in FirstStage_IV_SE]) + '\n')
f.write('Coef: ' + str([round(i, 3) for i in FirstStage_IV_Coef]
) + '\n')
f.write('t-stat: ' + str([round(i, 3) for i in FirstStage_IV_tstat]
) + '\n')
f.write('p-value: ' + str([round(i, 3) for i in FirstStage_IV_pvalue]
) + '\n')
VAR VARIABLE: IP cCPItaxSA JGB1Y TOPIXMA CBSm VAR SAMPLE: 1990-01-01 - 2020-01-01 IV VARIABLE: TARGET IV SAMPLE: 1992-07-01 - 2020-01-01 F-stat: 48.876 , robust F-stat: 15.11 R2: 0.129 Adj R2: 0.127 SE: [np.float64(0.004), np.float64(0.318)] Coef: [0.003, 1.235] t-stat: [0.668, 3.887] p-value: [0.504, 0.0]
（練習問題2）操作変数をEYF12とした場合¶
- 操作変数は、EYF12
- その他のセットアップは、本文と同様
In [ ]:
# データの読み込み
DATASET_VAR = pd.read_excel('/content/drive/My Drive/Data/data_VAR8.xlsx',
sheet_name='data_VAR', header=[0], index_col=0)
DATASET_IV = pd.read_excel('/content/drive/My Drive/Data/data_VAR8.xlsx',
sheet_name='data_IV', header=[0], index_col=0)
class structtype():
def __init__(self):
self.data = []
def psvar(VARresids, proxy, policyVarPos, nLag):
if np.sum(np.isnan(proxy)) > 0 and len(proxy) == len(VARresids):
index_notnan = ~np.isnan(proxy[:, 0])
proxy = proxy[index_notnan, :]
VARresids = VARresids[index_notnan, :]
VARresids = np.matrix(VARresids)
proxy = np.matrix(proxy)
T, N = VARresids.shape
t, n = proxy.shape
rows_iP = list(range(N))
rows_iP.remove(policyVarPos)
rows_iP = [policyVarPos] + rows_iP
rfResids = VARresids[:, rows_iP]
IV1X = np.concatenate([np.ones([t, 1]), proxy], axis=1)
IV1_coef = np.linalg.lstsq(IV1X, rfResids[:, 0], rcond=None)[0]
IV1_resids = rfResids[:, 0] - IV1X @ IV1_coef
IV_uphat = IV1X @ IV1_coef
IV2X = np.concatenate([np.ones([t, 1]), IV_uphat], axis=1)
IV2_coef = np.linalg.lstsq(IV2X, rfResids[:, 1:], rcond=None)[0]
IV2_resids = rfResids[:, 1:] - IV2X @ IV2_coef
s21is11 = IV2_coef[1:, :].T
SigmaU = np.cov(rfResids.T)
sigma11 = SigmaU[0:1, 0:1]
sigma12 = SigmaU[0:1, 1:]
sigma21 = SigmaU[1:, 0:1]
sigma22 = SigmaU[1:, 1:]
Qmat = np.dot(s21is11*sigma11, s21is11.T) - (
np.dot(sigma21, s21is11.T) + np.dot(s21is11, sigma21.T)) + sigma22
Qinv = np.linalg.solve(Qmat, (sigma21 - s21is11 * sigma11))
s12s12p = np.dot((sigma21- s21is11*sigma11).T, Qinv)
s11s11p = sigma11 - s12s12p
s11 = math.sqrt(s11s11p)
s1 = np.empty([N, 1])
s1[0, 0] = s11
s1[1:, 0:1] = s21is11*s11
IV1_dev = rfResids[:, 0] - np.mean(rfResids[:, 0])
SST_m = IV1_dev.T @ IV1_dev
SSE_m = IV1_resids[:, 0].T @ IV1_resids[:, 0]
Fstat = ((SST_m - SSE_m) / n) / (SSE_m / (t - n - 1))
R2 = 1 - SSE_m / SST_m
R2adj = R2 - (n / (t - n - 1)) * (1 - R2)
SS_m = np.zeros([n + 1, n + 1])
for ii in range(t):
SS_m = SS_m + (1 / t * IV1X[ii, :].T @ IV1X[ii, :] * (
IV1_resids[ii, 0] ** 2))
Avarb_m = np.linalg.inv(1 / t * IV1X.T @ IV1X) * SS_m * np.linalg.inv(
1 / t * IV1X.T @ IV1X)
IV1_robSE = np.sqrt(Avarb_m[1, 1])
IV1_robtstat = IV1_coef[1] / IV1_robSE
RR_m = np.concatenate([np.zeros([n, 1]), np.eye(n)], axis=1)
WW_m = t * (RR_m@IV1_coef[:, 0]).T @ np.linalg.inv(
RR_m @ Avarb_m @ RR_m.T) @ (RR_m @ IV1_coef[:, 0])
Fstat_rob = WW_m / n
res = structtype
res.B = s1[1:, :]
res.B = np.insert(res.B, policyVarPos, s1[0, 0], axis=0)
res.Coef = np.asarray(IV1_coef).reshape(-1).tolist()
res.Coef_const = np.asarray(IV1_coef[0, 0]).reshape(-1).tolist()
res.Obs = t
res.R2 = R2[0, 0]
res.R2adj = R2adj[0, 0]
res.Fstat = Fstat[0, 0]
res.Fstat_rob = Fstat_rob[0, 0]
Bzero = np.zeros([N, N - 1])
Bzero = np.insert(Bzero, policyVarPos, res.B.T, axis=1)
model_ols = OLS(rfResids[:,0], IV1X)
result = model_ols.fit(cov_type = 'HC0')
res.SE = result.bse
res.tstat = result.tvalues.tolist()
res.pvalue = result.pvalues.tolist()
result.fvalue
result.rsquared
result.rsquared_adj
return res, Bzero
VAR_varname = ['IP', 'cCPItaxSA', 'JGB1Y', 'TOPIXMA', 'CBSm']
VAR_figname = ['鉱工業生産', '消費者物価指数', '1年物国債利回り', '東証株価指数',
'社債スプレッド']
IV_varname = ['EYF12']
MPvar = 'JGB1Y'
VAR_smpl_from = '1990-01-01'
VAR_smpl_until = '2020-01-01'
IV_smpl_from = '1992-07-01'
IV_smpl_until = '2020-01-01'
iScheme = 'IV'
nLag = 12
nHor = 48
bandsCoverage = 0.9
shockSize = 1
def Gmatrices(AL, C, nLag, nHor, n):
J = np.hstack([I, np.zeros([n, (nLag-1)*n])])
Alut = np.vstack(
[AL, np.hstack([np.eye(n*(nLag-1)), np.zeros([n*(nLag-1), n])])])
AJ = np.zeros([n*nLag, n, nHor+1])
for j in range(nHor+1):
if j == 0:
AJ[:, :, j] = np.eye(Alut.shape[0]) @ J.T
else:
AJ[:, :, j] = np.linalg.matrix_power(Alut, j) @ J.T
JAp = np.reshape(AJ, [n*nLag, n*(nHor+1)], order='F').T
AJaux = np.zeros([JAp.shape[0]*n, JAp.shape[1]*n, nHor+1])
Caux = np.reshape(
np.hstack([I, C[:, :nHor*n]]), [n, n, nHor+1], order='F')
for i in range(nHor+1):
AJaux[((n**2)*i):, :, i] = np.kron(JAp[:n*(nHor+1-i), :],
Caux[:, :, i])
Gaux = np.transpose(np.reshape(
np.sum(AJaux, axis=2).T, [(n**2)*nLag, n**2, nHor+1], order='F'),
[1, 0, 2])
G = np.zeros([Gaux.shape[0], Gaux.shape[1], Gaux.shape[2]+1])
G[:, :, 1:] = Gaux
return Gaux
try:
VAR_smpl_from_index = DATASET_VAR.index.get_loc(VAR_smpl_from)
VAR_smpl_until_index = DATASET_VAR.index.get_loc(VAR_smpl_until)
except KeyError:
print('ERROR: CHECK THE VAR SAMPLE PERIOD YOU SET.')
if np.count_nonzero(np.isin(DATASET_VAR.columns, VAR_varname)) != len(
VAR_varname):
print('ERROR: CHECK THE VAR VARIABLE NAMES YOU SET.')
Ydf = DATASET_VAR.loc[VAR_smpl_from:VAR_smpl_until, VAR_varname]
print('VAR VARIABLE:', *Ydf.columns, sep=' ')
print('VAR SAMPLE: ', dt.date(Ydf.index[0]), '-', dt.date(Ydf.index[-1]))
if iScheme == 'IV':
IVstart = DATASET_IV[IV_varname].dropna().index[0]
if dt.strptime(IV_smpl_from, '%Y-%m-%d') < IVstart:
IV_smpl_from = dt.strftime(IVstart, '%Y-%m-%d')
try:
VAR_smpl_from_IV_index = DATASET_VAR.index.get_loc(IV_smpl_from)
VAR_smpl_until_IV_index = DATASET_VAR.index.get_loc(IV_smpl_until)
except KeyError:
print('ERROR: The IV sample is longer than the VAR sample period.')
if VAR_smpl_from_IV_index < VAR_smpl_from_index + nLag:
VAR_smpl_from_IV_index = VAR_smpl_from_index + nLag
IV_smpl_from = DATASET_VAR.index[VAR_smpl_from_IV_index]
if VAR_smpl_until_IV_index > VAR_smpl_until_index:
IV_smpl_until = VAR_smpl_until
VAR_smpl_until_IV_index = VAR_smpl_until_index
IVindex = np.empty([2, 1], dtype='int')
IVindex[0] = VAR_smpl_from_IV_index - (VAR_smpl_from_index + nLag)
IVindex[1] = VAR_smpl_until_index - VAR_smpl_until_IV_index
try:
IV_smpl_from_index = DATASET_IV.index.get_loc(IV_smpl_from)
IV_smpl_until_index = DATASET_IV.index.get_loc(IV_smpl_until)
except KeyError:
print('ERROR: CHECK THE IV SAMPLE PERIOD YOU SET.')
if np.count_nonzero(np.isin(DATASET_IV.columns, IV_varname)) != len(
IV_varname):
print('ERROR: CHECK THE INSTRUMENT NAME YOU SET.')
IVdf = DATASET_IV.loc[IV_smpl_from:IV_smpl_until, IV_varname]
print('IV VARIABLE:', *IVdf.columns, sep=' ')
print('IV SAMPLE: ', dt.date(IVdf.index[0]), '-', dt.date(IVdf.index[-1]))
elif iScheme == 'CHOL':
print('IDENTIFICATION: Cholesky decomposition')
if Ydf.iloc[0, :].isnull().sum() > 0:
print('*** ERROR: CHECK THE BEGINNING OF THE VAR SAMPLE PERIOD ***')
if Ydf.iloc[-1, :].isnull().sum() > 0:
print('*** ERROR: CHECK THE END OF THE VAR SAMPLE PERIOD ***')
Y = Ydf.values
if iScheme == 'IV':
IV = IVdf.values
MPshockIndex = VAR_varname.index(MPvar)
T, n = Y.shape
shockS = np.zeros([n, 1])
sLevel = abs(stats.norm.ppf((1-bandsCoverage)/2))
IRF = np.empty([n, nHor+1])
IRF.fill(np.nan)
Yh = Y[nLag:, :]
nT, nY = Yh.shape
YhLag = np.empty([T-nLag, n*nLag])
YhLag.fill(np.nan)
for j in range(nLag):
YhLag[:, n*j:n*(j+1)] = Y[nLag-j-1:-j-1, :].copy()
YhprojSet = np.concatenate([np.ones([nT, 1]), YhLag], axis=1)
projCoeffs = np.linalg.lstsq(YhprojSet, Yh, rcond=None)[0]
irfCoeffs = projCoeffs[1:n+1, :]
u = Yh - YhprojSet @ projCoeffs
SigmaU = np.cov(u.T)
if iScheme == 'CHOL':
Bzero = np.linalg.cholesky(SigmaU)
elif iScheme == 'IV':
if IVindex[1, 0] == 0:
FirstStage, Bzero = psvar(
u[IVindex[0,0]:, :], IV, MPshockIndex, nLag)
else:
FirstStage, Bzero = psvar(
u[IVindex[0, 0]:-IVindex[1, 0], :], IV, MPshockIndex, nLag)
print('F-stat: ', round(FirstStage.Fstat, 3),
', robust F-stat:', round(FirstStage.Fstat_rob, 3))
print('R2:', round(FirstStage.R2, 3))
print('Adj R2:', round(FirstStage.R2adj, 3))
print('SE:', [round(i, 3) for i in FirstStage.SE])
print('Coef:', [round(i, 3) for i in FirstStage.Coef])
print('t-stat:', [round(i, 3) for i in FirstStage.tstat])
print('p-value:', [round(i, 3) for i in FirstStage.pvalue])
FirstStage_IV_Fstat = FirstStage.Fstat.copy()
FirstStage_IV_Fstat_rob = FirstStage.Fstat_rob.copy()
FirstStage_IV_R2 = FirstStage.R2.copy()
FirstStage_IV_R2adj = FirstStage.R2adj.copy()
FirstStage_IV_SE = FirstStage.SE.copy()
FirstStage_IV_Coef = FirstStage.Coef.copy()
FirstStage_IV_tstat = FirstStage.tstat.copy()
FirstStage_IV_pvalue = FirstStage.pvalue.copy()
if shockSize == False:
shockS[MPshockIndex] = 1
else:
shockS[MPshockIndex] = shockSize / Bzero[MPshockIndex, MPshockIndex]
IRF[:, 0:1] = Bzero @ shockS
A = np.zeros([n*nLag, n*nLag])
A[0:n, :] = projCoeffs[1:, :].T
A[n:, 0:n*(nLag-1)] = np.eye(n*(nLag-1))
for h in range(nHor):
Ah = np.linalg.matrix_power(A, h+1)
IRF[:, h+1:h+2] = Ah[0:n, 0:n] @ Bzero @ shockS
m = YhprojSet.shape[1] - (n*nLag)
XSVARp = YhprojSet
if iScheme == 'IV':
k = IV.shape[1] + 1
IVm = np.empty([Y.shape[0], k])
IVm.fill(np.nan)
if IVindex[1, 0] == 0:
IVm[IVindex[0, 0]+nLag:, 0] = 1
IVm[IVindex[0, 0]+nLag:, 1:] = IV
else:
IVm[IVindex[0, 0]+nLag:-IVindex[1, 0], 0] = 1
IVm[IVindex[0, 0]+nLag:-IVindex[1, 0], 1:] = IV
IVm = IVm[nLag:, :]
IVm = np.nan_to_num(IVm)
else:
k=1
IVm = np.zeros([349, 1])
matagg = np.hstack([XSVARp, u, IVm]).T
T1aux = u.shape[0]
T2aux = matagg.shape[0]
etaaux = np.reshape(u.T, [nY, 1, T1aux], order='F')
mataggaux = np.transpose(np.reshape(matagg, [T2aux, 1, T1aux]), [1, 0, 2])
auxeta = np.multiply(etaaux, mataggaux) - \
np.mean(np.multiply(etaaux, mataggaux), 2)[:, :, np.newaxis]
vecAss1 = np.reshape(
auxeta, [(nY*m)+(nLag*(nY**2))+nY**2+(nY*k), 1, T1aux], order='F')
AuxHAC1 = vecAss1
AuxHAC2 = np.reshape(
AuxHAC1, [vecAss1.shape[0], vecAss1.shape[2]], order='F').T
NWlags = int(np.floor(4*((nT/100) ** (2/9))))
Sigma0 = (1/(AuxHAC2.shape[0]))*(AuxHAC2.T@AuxHAC2)
SigmaNW = Sigma0.copy()
for i in range(NWlags):
Sigma_cov_i = (1/AuxHAC2.shape[0]) * (AuxHAC2[:-i-1, :].T@AuxHAC2[1+i:, :])
SigmaNW = SigmaNW + (1-(i+1)/(NWlags+1)) * (Sigma_cov_i+Sigma_cov_i.T)
AuxHAC3 = SigmaNW
WhatAss1 = AuxHAC3
I = np.eye(nY)
V = np.kron(I[0, :], I)
for i_vars in range(1, n):
V = np.vstack([V, np.kron(I[i_vars, :], I[i_vars:, :])])
Q1 = XSVARp.T @ XSVARp / T1aux
Q2 = IVm.T @ XSVARp / T1aux
invQ1 = np.linalg.inv(Q1)
Shat_a1 = np.kron(np.hstack([np.zeros([n*nLag, m]), np.eye(n*nLag)])@invQ1, I)
Shat_a2 = np.zeros([(n**2)*nLag, n**2+(k*n)])
Shat_b = np.hstack([np.zeros([int(n*(n+1)/2), ((n**2)*nLag)+n *
m]), V, np.zeros([int(n*(n+1)/2), k*n])])
Shat_c = np.hstack(
[- np.kron(Q2@invQ1, I), np.zeros([k*n, n**2]), np.eye(k*n)])
Shat = np.vstack([np.hstack([Shat_a1, Shat_a2]), Shat_b, Shat_c])
WHataux = Shat @ WhatAss1 @ Shat.T
WHat11 = WHataux[:(n**2)*nLag, :(n**2)*nLag]
WHat12 = WHataux[:(n**2)*nLag, ((n**2)*nLag)+int(n*(n+1)/2):]
WHat22 = WHataux[((n**2)*nLag)+int(n*(n+1)/2):,
((n**2)*nLag)+int(n*(n+1)/2):]
WHatSigma = WHataux[(n**2)*nLag:((n**2)*nLag)+int(n*(n+1)/2),
(n**2)*nLag:((n**2)*nLag)+int(n*(n+1)/2)]
What = np.vstack([np.hstack([WHat11, WHat12]),
np.hstack([WHat12.T, WHat22])])
AL = projCoeffs[1:, :].T
vecAL = np.reshape(AL, [n, n, nLag], order='F')
vecALrevT = np.zeros([n, n, nHor])
for i in range(nHor):
if i < (nHor-nLag):
vecALrevT[:, :, i] = np.zeros([n, n])
else:
vecALrevT[:, :, i] = vecAL[:, :, nHor-i-1].T
vecALrevT = np.reshape(vecALrevT, [n, n*(nHor)], order='F')
MARepC = np.tile(vecAL[:, :, 0], [1, nHor])
for i in range(nHor-1):
MARepC[:, (n*(i+1)):(n*(i+2))] = np.hstack([I, MARepC[:, :n*(i+1)]]
) @ vecALrevT[:, (
nHor*n-(n*(i+2))):].T
Caux = np.hstack([I, MARepC])
C = np.reshape(Caux, [n, n, nHor+1], order='F')
G = Gmatrices(AL, Caux, nLag, nHor, n)
scale = 1
clevel = bandsCoverage
critval = sp.stats.norm.ppf(1-(1-clevel)/2, 0, 1) ** 2
lambdahat = np.zeros([n, nHor+1])
DmethodVar = np.zeros([n, nHor+1])
Dmethodlbound = np.zeros([n, nHor+1])
Dmethodubound = np.zeros([n, nHor+1])
if iScheme == 'IV':
Gamma_m = np.zeros([n, k])
for i in range(k):
Gamma_m[:,i] = u.T @ IVm[:,i] / nT
BzeroMP = Bzero[:, MPshockIndex][:, np.newaxis]
if shockSize == False:
Bzero00 = 1
else:
Bzero00 = shockSize*BzeroMP[MPshockIndex, 0]
for ih in range(nHor+1):
for ivar in range(n):
lambdahat[ivar, ih] = scale * \
I[:, ivar].T @ C[:, :, ih] @ BzeroMP / Bzero00
d1 = scale * np.kron(BzeroMP.T, I[:, ivar].T) @ G[:, :, ih]
sm = np.kron(Gamma_m[ivar,:], (I[:, ivar].T @ C[
:, :, ih] - lambdahat[ivar, ih] * I[:, MPshockIndex].T))
d = np.vstack([d1.T, sm[:, np.newaxis]])
DmethodVar[ivar, ih] = d.T @ What @ d
Dmethodlbound[ivar, ih] = lambdahat[ivar, ih] - \
((critval/nT)**.5)*(DmethodVar[ivar, ih]**.5)/abs(Bzero00)
Dmethodubound[ivar, ih] = lambdahat[ivar, ih] + \
((critval/nT)**.5)*(DmethodVar[ivar, ih]**.5)/abs(Bzero00)
else:
Bzero00 = Bzero[MPshockIndex, MPshockIndex]
w = np.arange(nY*nY).reshape((nY, nY), order='F').T.ravel(order='F')
Kkk = np.eye(nY*nY)[w, :]
H = V.T @ np.linalg.inv(V @ (np.eye(nY**2)+Kkk) @ np.kron(Bzero,I) @ V.T)
for ih in range(nHor+1):
lambdahat[:, ih] = scale * I @ C[:, :, ih] @ Bzero[
:, MPshockIndex] / Bzero00
lambdahat_all_ih = scale * I @ C[:, :, ih] @ Bzero
if ih == 0:
d1 = np.zeros([nY**2,nLag*nY**2])
else:
d1 = scale * np.kron(Bzero.T, I) @ G[:, :, ih]
d2 = scale * np.kron(I, lambdahat_all_ih) @ H
DmethodVar_h = d1 @ WHat11 @ d1.T + d2 @ WHatSigma @ d2.T
DmethodVar[:, ih] = np.diag(DmethodVar_h)[
(MPshockIndex-1)*nY:MPshockIndex*nY]
Dmethodlbound[:, ih] = lambdahat[:, ih] - \
((critval/nT)**.5)*(DmethodVar[:, ih]**.5) / Bzero00
Dmethodubound[:, ih] = lambdahat[:, ih] + \
((critval/nT)**.5)*(DmethodVar[:, ih]**.5) / Bzero00
IRFl = Dmethodlbound.copy()
IRFh = Dmethodubound.copy()
IRF = lambdahat
if iScheme == 'CHOL':
figcolor = 'r'
elif iScheme == 'IV':
figcolor = 'b'
if not 'VAR_figname' in locals():
VAR_figname = VAR_varname
fig_hrzn = np.arange(0, nHor+1)
nvars = len(VAR_varname)
fig = plt.figure(figsize=(nvars*3, 8))
ylim=([-10, 6], [-1.4, 0.4], [0, 1.8], [-50, 10], [-0.4, 0.6])
for i in range(nvars):
ax = fig.add_subplot(2, 3, i+1)
ax.set_title(VAR_figname[i], fontsize=14)
ax.axhline(y=0, xmin=0, xmax=1, color='black', linewidth=0.75)
ax.plot(fig_hrzn, IRF[i, :], color='black', linewidth=2)
ax.plot(fig_hrzn, IRFl[i, :], color='black', linewidth=1,
linestyle='dashed')
ax.plot(fig_hrzn, IRFh[i, :], color='black', linewidth=1,
linestyle='dashed')
ax.set_xlim([0, nHor])
ax.set_xlabel('月', fontsize=14)
ax.set_xticks(range(0, nHor+1, 6))
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.set_ylim(ylim[i])
fig.tight_layout()
SaveDirName = 'Figure/VAR'
if not os.path.exists('Figure'):
os.makedirs('Figure')
if not os.path.exists(SaveDirName):
os.makedirs(SaveDirName)
if iScheme == 'CHOL':
FigName = 'irsSVAR_JP_Chol'
elif iScheme == 'IV':
FigName = 'irsSVAR_JP_IV'
SaveFigName = SaveDirName + '/' + FigName
with open(SaveFigName+'.txt', 'w') as f:
f.write('LAG ORDER: ' + str(nLag) + '\n')
f.write('POLICY INDICATOR: ' + MPvar + '\n')
f.write('VAR VARIABLE: ' + str(Ydf.columns.values) + '\n')
f.write('VAR SAMPLE: ' + str(dt.date(Ydf.index[0])) + ' - ' + str(
dt.date(Ydf.index[-1])) + '\n')
if iScheme == 'IV':
f.write('IV VARIABLE: ' + str(IVdf.columns.values) + '\n')
f.write('IV SAMPLE: ' + str(dt.date(IVdf.index[0])) + ' - ' + str(
dt.date(IVdf.index[-1])) + '\n')
f.write('F-stat: ' + str(round(FirstStage_IV_Fstat, 3)
) + ', robust F-stat: ' + str(round(FirstStage_IV_Fstat_rob, 3)) + '\n')
f.write('R2: ' + str(round(FirstStage_IV_R2, 3)) + '\n')
f.write('Adj R2: ' + str(round(FirstStage_IV_R2adj, 3)) + '\n')
f.write('SE: ' + str([round(i, 3) for i in FirstStage_IV_SE]) + '\n')
f.write('Coef: ' + str([round(i, 3) for i in FirstStage_IV_Coef]
) + '\n')
f.write('t-stat: ' + str([round(i, 3) for i in FirstStage_IV_tstat]
) + '\n')
f.write('p-value: ' + str([round(i, 3) for i in FirstStage_IV_pvalue]
) + '\n')
VAR VARIABLE: IP cCPItaxSA JGB1Y TOPIXMA CBSm VAR SAMPLE: 1990-01-01 - 2020-01-01 IV VARIABLE: EYF12 IV SAMPLE: 1992-07-01 - 2020-01-01 F-stat: 33.372 , robust F-stat: 11.519 R2: 0.092 Adj R2: 0.089 SE: [np.float64(0.004), np.float64(0.194)] Coef: [0.004, 0.66] t-stat: [1.014, 3.394] p-value: [0.311, 0.001]
