Scarica Soluzione esercizio econometria Python: calcolo statistiche regressione lineare e più Esercizi in PDF di Econometria solo su Docsity!
import numpy as np import statsmodels.api as sm from matplotlib import pyplot import matplotlib.pyplot as plt import math import pandas as pd import os import datetime import csv from scipy.ndimage.interpolation import shift import random from datetime import datetime from scipy.stats import norm from scipy.stats import f from scipy.stats import chi random.seed(datetime.now()) df1 = csv.reader(open('canada_cons_inc.txt'),delimiter='\t') b1 = np.array(list(df1))[ 1 :, 0 : 2 ].astype('float') nrow = b1.shape[ 0 ] n=nrow - 1 intercept = np.ones( (nrow- 1 , 1 ) ) C = np.zeros(n) C=b1[ 1 :, 0 ].reshape(- 1 , 1 ) Y = np.zeros(n) Y=b1[ 1 :, 1 ].reshape(- 1 , 1 ) c_lag = np.zeros(n) c_lag = b1[ 0 : 113 , 0 ].reshape(- 1 , 1 ) c_lag = c_lag[ 0 : 112 ] y_lag = np.zeros(n) y_lag = b1[ 0 : 113 , 1 ].reshape(- 1 , 1 ) y_lag = y_lag[ 0 : 112 ] x = np.concatenate((intercept,c_lag,Y,y_lag), axis= 1 ) beta = np.dot(np.linalg.inv(np.dot(x.T,x)),np.dot(x.T,C)) #print(beta) p=np.dot(np.dot(x,np.linalg.inv(np.dot(x.T,x))),x.T) c_fit=np.dot(np.dot(np.dot(x,np.linalg.inv(np.dot(x.T,x))),x.T),C) ident = np.identity(n) m=ident-p res = np.dot(m,C) asigma = np.dot(np.dot(C.T,m),C)/ 108 avar = asigma*np.linalg.inv(np.dot(x.T,x)) R = [ 0 , 0 , 1 , 1 ] R = np.array(R) RavarR = np.dot(np.dot(R,avar),R.T) RavarR_inv = 1 /RavarR Rb = np.dot(R,beta) wald = np.dot(np.dot(Rb.T,RavarR_inv),Rb) t_stat= np.sqrt(wald) #t-stat computed as Wald/F print("The t-stat is ",t_stat) print("The asymp pvalue is ", 2 ( 1 - norm.cdf(t_stat, 0 , 1 )))#Asymp. p-value print("The F-stat is ", wald) print("The pvalue is ", 1 - f.cdf(wald, 1 , 108 )) asigma = np.dot(np.dot(C.T,m),C)/ 112 avar = asigmanp.linalg.inv(np.dot(x.T,x)) R = [ 0 , 0 , 1 , 1 ]
R = np.array(R) RavarR = np.dot(np.dot(R,avar),R.T) RavarR_inv = 1 /RavarR Rb = np.dot(R,beta) wald = np.dot(np.dot(Rb.T,RavarR_inv),Rb) print("The asymp wald is ",wald) print("The pvalue is ", 1 - chi2.cdf(wald, 1 )) #print(wald) #EconFin_ex5rev20sol #********************************* #R = [0,0,1,0] #R = np.array(R) #RavarR = np.dot(np.dot(R,avar),R.T) #RavarR_inv = 1/RavarR #print(RavarR_inv) #Rb = np.dot(R,beta) #wald = np.dot(np.dot(Rb.T,RavarR_inv),Rb) #print(wald) #********************************* nx = y_lag - Y x_rest = np.concatenate((intercept,c_lag,nx), axis= 1 ) beta_r = np.dot(np.linalg.inv(np.dot(x_rest.T,x_rest)),np.dot(x_rest.T,C)) #print(beta_r) p_rest=np.dot(np.dot(x_rest,np.linalg.inv(np.dot(x_rest.T,x_rest))),x_rest.T) c_fit_rest=np.dot(np.dot(np.dot(x_rest,np.linalg.inv(np.dot(x_rest.T,x_rest))),x_rest.T),C) ident = np.identity(n) m_rest=ident-p_rest res_r = np.dot(m_rest,C) asigma_r = np.dot(np.dot(C.T,m_rest),C)/ 109 #print(res_r) y_b=np.zeros(n) est_std=np.sqrt(asigma_r) count_b= 0 b_reps= 999 for b in range( 0 ,b_reps): for i in range( 0 , 112 ): b_res_r=np.random.normal( 0 ,est_std) y_b[i]=beta_r[ 0 ]+beta_r[ 1 ]c_lag[i]+beta_r[ 2 ]nx[i]+b_res_r beta = np.dot(np.linalg.inv(np.dot(x.T,x)),np.dot(x.T,y_b)) p=np.dot(np.dot(x,np.linalg.inv(np.dot(x.T,x))),x.T) c_fit=np.dot(np.dot(np.dot(x,np.linalg.inv(np.dot(x.T,x))),x.T),y_b) ident = np.identity(n) m=ident-p res = np.dot(m,y_b) asigma = np.dot(np.dot(y_b.T,m),y_b)/ 108 avar = asigma*np.linalg.inv(np.dot(x.T,x)) R = [ 0 , 0 , 1 , 1 ] R = np.array(R) RavarR = np.dot(np.dot(R,avar),R.T) RavarR_inv = 1 /RavarR Rb = np.dot(R,beta) b_wald = np.dot(np.dot(Rb.T,RavarR_inv),Rb) b_t_stat=np.sqrt(b_wald) if b_t_stat > t_stat: count_b = count_b + 1 print("The parametric bootstrap p-vale is ", count_b/b_reps)
