[docs]defpredict(self,params,start=None,end=None,lags=1,trend='c'):""" Returns in-sample predictions or forecasts """start=self._get_predict_start(start,lags)end,out_of_sample=self._get_predict_end(end)ifend<start:raiseValueError("end is before start")ifend==start+out_of_sample:returnnp.array([])k_trend=util.get_trendorder(trend)k=self.neqsk_ar=lagspredictedvalues=np.zeros((end+1-start+out_of_sample,k))ifk_trend!=0:intercept=params[:k_trend]predictedvalues+=intercepty=self.yX=util.get_var_endog(y,lags,trend=trend,has_constant='raise')fittedvalues=np.dot(X,params)fv_start=start-k_arpv_end=min(len(predictedvalues),len(fittedvalues)-fv_start)fv_end=min(len(fittedvalues),end-k_ar+1)predictedvalues[:pv_end]=fittedvalues[fv_start:fv_end]ifnotout_of_sample:returnpredictedvalues# fit out of sampley=y[-k_ar:]coefs=params[k_trend:].reshape((k_ar,k,k)).swapaxes(1,2)predictedvalues[pv_end:]=forecast(y,coefs,intercept,out_of_sample)returnpredictedvalues

[docs]defplot(self):"""Plot input time series """plotting.plot_mts(self.y,names=self.names,index=self.dates)

@propertydefdf_model(self):"""Number of estimated parameters, including the intercept / trends """returnself.neqs*self.k_ar+self.k_trend@propertydefdf_resid(self):"""Number of observations minus number of estimated parameters"""returnself.nobs-self.df_model@cache_readonly

[docs]deffittedvalues(self):"""The predicted insample values of the response variables of the model. """returnnp.dot(self.ys_lagged,self.params)

def_omega_forc_cov(self,steps):# Approximate MSE matrix \Omega(h) as defined in Lut p97G=self._zzGinv=L.inv(G)# memoize powers of B for speedup# TODO: see if can memoize betterB=self._bmat_forc_cov()_B={}defbpow(i):ifinotin_B:_B[i]=np.linalg.matrix_power(B,i)return_B[i]phis=self.ma_rep(steps)sig_u=self.sigma_uomegas=np.zeros((steps,self.neqs,self.neqs))forhinrange(1,steps+1):ifh==1:omegas[h-1]=self.df_model*self.sigma_ucontinueom=omegas[h-1]foriinrange(h):forjinrange(h):Bi=bpow(h-1-i)Bj=bpow(h-1-j)mult=np.trace(chain_dot(Bi.T,Ginv,Bj,G))om+=mult*chain_dot(phis[i],sig_u,phis[j].T)omegas[h-1]=omreturnomegasdef_bmat_forc_cov(self):# B as defined on p. 96 of Lutupper=np.zeros((1,self.df_model))upper[0,0]=1lower_dim=self.neqs*(self.k_ar-1)I=np.eye(lower_dim)lower=np.column_stack((np.zeros((lower_dim,1)),I,np.zeros((lower_dim,self.neqs))))returnnp.vstack((upper,self.params.T,lower))

[docs]defreorder(self,order):"""Reorder variables for structural specification """iflen(order)!=len(self.params[0,:]):raiseValueError("Reorder specification length should match number of endogenous variables")#This convert order to list of integers if given as stringsifisinstance(order[0],string_types):order_new=[]fori,naminenumerate(order):order_new.append(self.names.index(order[i]))order=order_newreturn_reordered(self,order)

[docs]deftest_causality(self,equation,variables,kind='f',signif=0.05,verbose=True):"""Compute test statistic for null hypothesis of Granger-noncausality, general function to test joint Granger-causality of multiple variables Parameters ---------- equation : string or int Equation to test for causality variables : sequence (of strings or ints) List, tuple, etc. of variables to test for Granger-causality kind : {'f', 'wald'} Perform F-test or Wald (chi-sq) test signif : float, default 5% Significance level for computing critical values for test, defaulting to standard 0.95 level Notes ----- Null hypothesis is that there is no Granger-causality for the indicated variables. The degrees of freedom in the F-test are based on the number of variables in the VAR system, that is, degrees of freedom are equal to the number of equations in the VAR times degree of freedom of a single equation. Returns ------- results : dict """ifisinstance(variables,(string_types,int,np.integer)):variables=[variables]k,p=self.neqs,self.k_ar# number of restrictionsN=len(variables)*self.k_ar# Make restriction matrixC=np.zeros((N,k**2*p+k),dtype=float)eq_index=self.get_eq_index(equation)vinds=mat([self.get_eq_index(v)forvinvariables])# remember, vec is column order!offsets=np.concatenate([k+k**2*j+k*vinds+eq_indexforjinrange(p)])C[np.arange(N),offsets]=1# Lutkepohl 3.6.5Cb=np.dot(C,vec(self.params.T))middle=L.inv(chain_dot(C,self.cov_params,C.T))# wald statisticlam_wald=statistic=chain_dot(Cb,middle,Cb)ifkind.lower()=='wald':df=Ndist=stats.chi2(df)elifkind.lower()=='f':statistic=lam_wald/Ndf=(N,k*self.df_resid)dist=stats.f(*df)else:raiseException('kind %s not recognized'%kind)pvalue=dist.sf(statistic)crit_value=dist.ppf(1-signif)conclusion='fail to reject'ifstatistic<crit_valueelse'reject'results={'statistic':statistic,'crit_value':crit_value,'pvalue':pvalue,'df':df,'conclusion':conclusion,'signif':signif}ifverbose:summ=output.causality_summary(results,variables,equation,kind)print(summ)returnresults

classVARResultsWrapper(wrap.ResultsWrapper):_attrs={'bse':'columns_eq','cov_params':'cov','params':'columns_eq','pvalues':'columns_eq','tvalues':'columns_eq','sigma_u':'cov_eq','sigma_u_mle':'cov_eq','stderr':'columns_eq'}_wrap_attrs=wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_attrs,_attrs)_methods={}_wrap_methods=wrap.union_dicts(tsbase.TimeSeriesResultsWrapper._wrap_methods,_methods)_wrap_methods.pop('cov_params')# not yet a method in VARResultswrap.populate_wrapper(VARResultsWrapper,VARResults)

[docs]defsummary(self):buf=StringIO()rng=lrange(self.periods)foriinrange(self.neqs):ppm=output.pprint_matrix(self.decomp[i],rng,self.names)buf.write('FEVD for %s\n'%self.names[i])buf.write(ppm+'\n')print(buf.getvalue())

[docs]defplot(self,periods=None,figsize=(10,10),**plot_kwds):"""Plot graphical display of FEVD Parameters ---------- periods : int, default None Defaults to number originally specified. Can be at most that number """importmatplotlib.pyplotaspltk=self.neqsperiods=periodsorself.periodsfig,axes=plt.subplots(nrows=k,figsize=figsize)fig.suptitle('Forecast error variance decomposition (FEVD)')colors=[str(c)forcinnp.arange(k,dtype=float)/k]ticks=np.arange(periods)limits=self.decomp.cumsum(2)foriinrange(k):ax=axes[i]this_limits=limits[i].Thandles=[]forjinrange(k):lower=this_limits[j-1]ifj>0else0upper=this_limits[j]handle=ax.bar(ticks,upper-lower,bottom=lower,color=colors[j],label=self.names[j],**plot_kwds)handles.append(handle)ax.set_title(self.names[i])# just use the last axis to get handles for plottinghandles,labels=ax.get_legend_handles_labels()fig.legend(handles,labels,loc='upper right')plotting.adjust_subplots(right=0.85)