# -*- coding: utf-8 -*-
"""
Impulse reponse-related code
"""
from __future__ import division
import numpy as np
import numpy.linalg as la
import scipy.linalg as L
from scipy import stats
from statsmodels.tools.decorators import cache_readonly
from statsmodels.tools.tools import chain_dot
#from statsmodels.tsa.api import VAR
from statsmodels.compat.python import range
import statsmodels.tsa.tsatools as tsa
import statsmodels.tsa.vector_ar.plotting as plotting
import statsmodels.tsa.vector_ar.util as util
mat = np.array
class BaseIRAnalysis(object):
"""
Base class for plotting and computing IRF-related statistics, want to be
able to handle known and estimated processes
"""
def __init__(self, model, P=None, periods=10, order=None, svar=False,
vecm=False):
self.model = model
self.periods = periods
self.neqs, self.lags, self.T = model.neqs, model.k_ar, model.nobs
self.order = order
if P is None:
sigma = model.sigma_u
# TODO, may be difficult at the moment
# if order is not None:
# indexer = [model.get_eq_index(name) for name in order]
# sigma = sigma[:, indexer][indexer, :]
# if sigma.shape != model.sigma_u.shape:
# raise ValueError('variable order is wrong length')
P = la.cholesky(sigma)
self.P = P
self.svar = svar
self.irfs = model.ma_rep(periods)
if svar:
self.svar_irfs = model.svar_ma_rep(periods, P=P)
else:
self.orth_irfs = model.orth_ma_rep(periods, P=P)
self.cum_effects = self.irfs.cumsum(axis=0)
if svar:
self.svar_cum_effects = self.svar_irfs.cumsum(axis=0)
else:
self.orth_cum_effects = self.orth_irfs.cumsum(axis=0)
# long-run effects may be infinite for VECMs.
if not vecm:
self.lr_effects = model.long_run_effects()
if svar:
self.svar_lr_effects = np.dot(model.long_run_effects(), P)
else:
self.orth_lr_effects = np.dot(model.long_run_effects(), P)
# auxiliary stuff
if vecm:
self._A = util.comp_matrix(model.var_rep)
else:
self._A = util.comp_matrix(model.coefs)
def cov(self, *args, **kwargs):
raise NotImplementedError
def cum_effect_cov(self, *args, **kwargs):
raise NotImplementedError
def plot(self, orth=False, impulse=None, response=None,
signif=0.05, plot_params=None, subplot_params=None,
plot_stderr=True, stderr_type='asym', repl=1000,
seed=None, component=None):
"""
Plot impulse responses
Parameters
----------
orth : bool, default False
Compute orthogonalized impulse responses
impulse : string or int
variable providing the impulse
response : string or int
variable affected by the impulse
signif : float (0 < signif < 1)
Significance level for error bars, defaults to 95% CI
subplot_params : dict
To pass to subplot plotting funcions. Example: if fonts are too big,
pass {'fontsize' : 8} or some number to your taste.
plot_params : dict
plot_stderr: bool, default True
Plot standard impulse response error bands
stderr_type: string
'asym': default, computes asymptotic standard errors
'mc': monte carlo standard errors (use rpl)
repl: int, default 1000
Number of replications for Monte Carlo and Sims-Zha standard errors
seed: int
np.random.seed for Monte Carlo replications
component: array or vector of principal component indices
"""
periods = self.periods
model = self.model
svar = self.svar
if orth and svar:
raise ValueError("For SVAR system, set orth=False")
if orth:
title = 'Impulse responses (orthogonalized)'
irfs = self.orth_irfs
elif svar:
title = 'Impulse responses (structural)'
irfs = self.svar_irfs
else:
title = 'Impulse responses'
irfs = self.irfs
if plot_stderr == False:
stderr = None
elif stderr_type not in ['asym', 'mc', 'sz1', 'sz2','sz3']:
raise ValueError("Error type must be either 'asym', 'mc','sz1','sz2', or 'sz3'")
else:
if stderr_type == 'asym':
stderr = self.cov(orth=orth)
if stderr_type == 'mc':
stderr = self.errband_mc(orth=orth, svar=svar,
repl=repl, signif=signif,
seed=seed)
if stderr_type == 'sz1':
stderr = self.err_band_sz1(orth=orth, svar=svar,
repl=repl, signif=signif,
seed=seed,
component=component)
if stderr_type == 'sz2':
stderr = self.err_band_sz2(orth=orth, svar=svar,
repl=repl, signif=signif,
seed=seed,
component=component)
if stderr_type == 'sz3':
stderr = self.err_band_sz3(orth=orth, svar=svar,
repl=repl, signif=signif,
seed=seed,
component=component)
fig = plotting.irf_grid_plot(irfs, stderr, impulse, response,
self.model.names, title, signif=signif,
subplot_params=subplot_params,
plot_params=plot_params,
stderr_type=stderr_type)
return fig
def plot_cum_effects(self, orth=False, impulse=None, response=None,
signif=0.05, plot_params=None,
subplot_params=None, plot_stderr=True,
stderr_type='asym', repl=1000, seed=None):
"""
Plot cumulative impulse response functions
Parameters
----------
orth : bool, default False
Compute orthogonalized impulse responses
impulse : string or int
variable providing the impulse
response : string or int
variable affected by the impulse
signif : float (0 < signif < 1)
Significance level for error bars, defaults to 95% CI
subplot_params : dict
To pass to subplot plotting funcions. Example: if fonts are too big,
pass {'fontsize' : 8} or some number to your taste.
plot_params : dict
plot_stderr: bool, default True
Plot standard impulse response error bands
stderr_type: string
'asym': default, computes asymptotic standard errors
'mc': monte carlo standard errors (use rpl)
repl: int, default 1000
Number of replications for monte carlo standard errors
seed: int
np.random.seed for Monte Carlo replications
"""
if orth:
title = 'Cumulative responses responses (orthogonalized)'
cum_effects = self.orth_cum_effects
lr_effects = self.orth_lr_effects
else:
title = 'Cumulative responses'
cum_effects = self.cum_effects
lr_effects = self.lr_effects
if stderr_type not in ['asym', 'mc']:
raise TypeError
else:
if stderr_type == 'asym':
stderr = self.cum_effect_cov(orth=orth)
if stderr_type == 'mc':
stderr = self.cum_errband_mc(orth=orth, repl=repl,
signif=signif, seed=seed)
if not plot_stderr:
stderr = None
fig = plotting.irf_grid_plot(cum_effects, stderr, impulse, response,
self.model.names, title, signif=signif,
hlines=lr_effects,
subplot_params=subplot_params,
plot_params=plot_params,
stderr_type=stderr_type)
return fig
[docs]class IRAnalysis(BaseIRAnalysis):
"""
Impulse response analysis class. Computes impulse responses, asymptotic
standard errors, and produces relevant plots
Parameters
----------
model : VAR instance
Notes
-----
Using Lütkepohl (2005) notation
"""
def __init__(self, model, P=None, periods=10, order=None, svar=False,
vecm=False):
BaseIRAnalysis.__init__(self, model, P=P, periods=periods,
order=order, svar=svar, vecm=vecm)
if vecm:
self.cov_a = model.cov_var_repr
else:
self.cov_a = model._cov_alpha
self.cov_sig = model._cov_sigma
# memoize dict for G matrix function
self._g_memo = {}
[docs] def cov(self, orth=False):
"""
Compute asymptotic standard errors for impulse response coefficients
Notes
-----
Lütkepohl eq 3.7.5
Returns
-------
"""
if orth:
return self._orth_cov()
covs = self._empty_covm(self.periods + 1)
covs[0] = np.zeros((self.neqs ** 2, self.neqs ** 2))
for i in range(1, self.periods + 1):
Gi = self.G[i - 1]
covs[i] = chain_dot(Gi, self.cov_a, Gi.T)
return covs
[docs] def errband_mc(self, orth=False, svar=False, repl=1000,
signif=0.05, seed=None, burn=100):
"""
IRF Monte Carlo integrated error bands
"""
model = self.model
periods = self.periods
if svar == True:
return model.sirf_errband_mc(orth=orth, repl=repl, T=periods,
signif=signif, seed=seed,
burn=burn, cum=False)
else:
return model.irf_errband_mc(orth=orth, repl=repl, T=periods,
signif=signif, seed=seed,
burn=burn, cum=False)
[docs] def err_band_sz1(self, orth=False, svar=False, repl=1000,
signif=0.05, seed=None, burn=100, component=None):
"""
IRF Sims-Zha error band method 1. Assumes symmetric error bands around
mean.
Parameters
----------
orth : bool, default False
Compute orthogonalized impulse responses
repl : int, default 1000
Number of MC replications
signif : float (0 < signif < 1)
Significance level for error bars, defaults to 95% CI
seed : int, default None
np.random seed
burn : int, default 100
Number of initial simulated obs to discard
component : neqs x neqs array, default to largest for each
Index of column of eigenvector/value to use for each error band
Note: period of impulse (t=0) is not included when computing
principle component
References
----------
Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse
Response". Econometrica 67: 1113-1155.
"""
model = self.model
periods = self.periods
if orth:
irfs = self.orth_irfs
elif svar:
irfs = self.svar_irfs
else:
irfs = self.irfs
neqs = self.neqs
irf_resim = model.irf_resim(orth=orth, repl=repl, T=periods, seed=seed,
burn=100)
q = util.norm_signif_level(signif)
W, eigva, k =self._eigval_decomp_SZ(irf_resim)
if component != None:
if np.shape(component) != (neqs,neqs):
raise ValueError("Component array must be " + str(neqs) + " x " + str(neqs))
if np.argmax(component) >= neqs*periods:
raise ValueError("Atleast one of the components does not exist")
else:
k = component
# here take the kth column of W, which we determine by finding the largest eigenvalue of the covaraince matrix
lower = np.copy(irfs)
upper = np.copy(irfs)
for i in range(neqs):
for j in range(neqs):
lower[1:,i,j] = irfs[1:,i,j] + W[i,j,:,k[i,j]]*q*np.sqrt(eigva[i,j,k[i,j]])
upper[1:,i,j] = irfs[1:,i,j] - W[i,j,:,k[i,j]]*q*np.sqrt(eigva[i,j,k[i,j]])
return lower, upper
[docs] def err_band_sz2(self, orth=False, svar=False, repl=1000, signif=0.05,
seed=None, burn=100, component=None):
"""
IRF Sims-Zha error band method 2.
This method Does not assume symmetric error bands around mean.
Parameters
----------
orth : bool, default False
Compute orthogonalized impulse responses
repl : int, default 1000
Number of MC replications
signif : float (0 < signif < 1)
Significance level for error bars, defaults to 95% CI
seed : int, default None
np.random seed
burn : int, default 100
Number of initial simulated obs to discard
component : neqs x neqs array, default to largest for each
Index of column of eigenvector/value to use for each error band
Note: period of impulse (t=0) is not included when computing
principle component
References
----------
Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse
Response". Econometrica 67: 1113-1155.
"""
model = self.model
periods = self.periods
if orth:
irfs = self.orth_irfs
elif svar:
irfs = self.svar_irfs
else:
irfs = self.irfs
neqs = self.neqs
irf_resim = model.irf_resim(orth=orth, repl=repl, T=periods, seed=seed,
burn=100)
W, eigva, k = self._eigval_decomp_SZ(irf_resim)
if component != None:
if np.shape(component) != (neqs,neqs):
raise ValueError("Component array must be " + str(neqs) + " x " + str(neqs))
if np.argmax(component) >= neqs*periods:
raise ValueError("Atleast one of the components does not exist")
else:
k = component
gamma = np.zeros((repl, periods+1, neqs, neqs))
for p in range(repl):
for i in range(neqs):
for j in range(neqs):
gamma[p,1:,i,j] = W[i,j,k[i,j],:] * irf_resim[p,1:,i,j]
gamma_sort = np.sort(gamma, axis=0) #sort to get quantiles
indx = round(signif/2*repl)-1,round((1-signif/2)*repl)-1
lower = np.copy(irfs)
upper = np.copy(irfs)
for i in range(neqs):
for j in range(neqs):
lower[:,i,j] = irfs[:,i,j] + gamma_sort[indx[0],:,i,j]
upper[:,i,j] = irfs[:,i,j] + gamma_sort[indx[1],:,i,j]
return lower, upper
[docs] def err_band_sz3(self, orth=False, svar=False, repl=1000, signif=0.05,
seed=None, burn=100, component=None):
"""
IRF Sims-Zha error band method 3. Does not assume symmetric error bands around mean.
Parameters
----------
orth : bool, default False
Compute orthogonalized impulse responses
repl : int, default 1000
Number of MC replications
signif : float (0 < signif < 1)
Significance level for error bars, defaults to 95% CI
seed : int, default None
np.random seed
burn : int, default 100
Number of initial simulated obs to discard
component : vector length neqs, default to largest for each
Index of column of eigenvector/value to use for each error band
Note: period of impulse (t=0) is not included when computing
principle component
References
----------
Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse
Response". Econometrica 67: 1113-1155.
"""
model = self.model
periods = self.periods
if orth:
irfs = self.orth_irfs
elif svar:
irfs = self.svar_irfs
else:
irfs = self.irfs
neqs = self.neqs
irf_resim = model.irf_resim(orth=orth, repl=repl, T=periods, seed=seed,
burn=100)
stack = np.zeros((neqs, repl, periods*neqs))
#stack left to right, up and down
for p in range(repl):
for i in range(neqs):
stack[i, p,:] = np.ravel(irf_resim[p,1:,:,i].T)
stack_cov=np.zeros((neqs, periods*neqs, periods*neqs))
W = np.zeros((neqs, periods*neqs, periods*neqs))
eigva = np.zeros((neqs, periods*neqs))
k = np.zeros((neqs))
if component != None:
if np.size(component) != (neqs):
raise ValueError("Component array must be of length " + str(neqs))
if np.argmax(component) >= neqs*periods:
raise ValueError("Atleast one of the components does not exist")
else:
k = component
#compute for eigen decomp for each stack
for i in range(neqs):
stack_cov[i] = np.cov(stack[i],rowvar=0)
W[i], eigva[i], k[i] = util.eigval_decomp(stack_cov[i])
gamma = np.zeros((repl, periods+1, neqs, neqs))
for p in range(repl):
c=0
for j in range(neqs):
for i in range(neqs):
gamma[p,1:,i,j] = W[j,k[j],i*periods:(i+1)*periods] * irf_resim[p,1:,i,j]
if i == neqs-1:
gamma[p,1:,i,j] = W[j,k[j],i*periods:] * irf_resim[p,1:,i,j]
gamma_sort = np.sort(gamma, axis=0) #sort to get quantiles
indx = round(signif/2*repl)-1,round((1-signif/2)*repl)-1
lower = np.copy(irfs)
upper = np.copy(irfs)
for i in range(neqs):
for j in range(neqs):
lower[:,i,j] = irfs[:,i,j] + gamma_sort[indx[0],:,i,j]
upper[:,i,j] = irfs[:,i,j] + gamma_sort[indx[1],:,i,j]
return lower, upper
def _eigval_decomp_SZ(self, irf_resim):
"""
Returns
-------
W: array of eigenvectors
eigva: list of eigenvalues
k: matrix indicating column # of largest eigenvalue for each c_i,j
"""
neqs = self.neqs
periods = self.periods
cov_hold = np.zeros((neqs, neqs, periods, periods))
for i in range(neqs):
for j in range(neqs):
cov_hold[i,j,:,:] = np.cov(irf_resim[:,1:,i,j],rowvar=0)
W = np.zeros((neqs, neqs, periods, periods))
eigva = np.zeros((neqs, neqs, periods, 1))
k = np.zeros((neqs, neqs))
for i in range(neqs):
for j in range(neqs):
W[i,j,:,:], eigva[i,j,:,0], k[i,j] = util.eigval_decomp(cov_hold[i,j,:,:])
return W, eigva, k
[docs] @cache_readonly
def G(self):
# Gi matrices as defined on p. 111
K = self.neqs
# nlags = self.model.p
# J = np.hstack((np.eye(K),) + (np.zeros((K, K)),) * (nlags - 1))
def _make_g(i):
# p. 111 Lutkepohl
G = 0.
for m in range(i):
# be a bit cute to go faster
idx = i - 1 - m
if idx in self._g_memo:
apow = self._g_memo[idx]
else:
apow = la.matrix_power(self._A.T, idx)
# apow = np.dot(J, apow)
apow = apow[:K]
self._g_memo[idx] = apow
# take first K rows
piece = np.kron(apow, self.irfs[m])
G = G + piece
return G
return [_make_g(i) for i in range(1, self.periods + 1)]
def _orth_cov(self):
# Lutkepohl 3.7.8
Ik = np.eye(self.neqs)
PIk = np.kron(self.P.T, Ik)
H = self.H
covs = self._empty_covm(self.periods + 1)
for i in range(self.periods + 1):
if i == 0:
apiece = 0
else:
Ci = np.dot(PIk, self.G[i-1])
apiece = chain_dot(Ci, self.cov_a, Ci.T)
Cibar = np.dot(np.kron(Ik, self.irfs[i]), H)
bpiece = chain_dot(Cibar, self.cov_sig, Cibar.T) / self.T
# Lutkepohl typo, cov_sig correct
covs[i] = apiece + bpiece
return covs
[docs] def cum_effect_cov(self, orth=False):
"""
Compute asymptotic standard errors for cumulative impulse response
coefficients
Parameters
----------
orth : boolean
Notes
-----
eq. 3.7.7 (non-orth), 3.7.10 (orth)
Returns
-------
"""
Ik = np.eye(self.neqs)
PIk = np.kron(self.P.T, Ik)
F = 0.
covs = self._empty_covm(self.periods + 1)
for i in range(self.periods + 1):
if i > 0:
F = F + self.G[i - 1]
if orth:
if i == 0:
apiece = 0
else:
Bn = np.dot(PIk, F)
apiece = chain_dot(Bn, self.cov_a, Bn.T)
Bnbar = np.dot(np.kron(Ik, self.cum_effects[i]), self.H)
bpiece = chain_dot(Bnbar, self.cov_sig, Bnbar.T) / self.T
covs[i] = apiece + bpiece
else:
if i == 0:
covs[i] = np.zeros((self.neqs**2, self.neqs**2))
continue
covs[i] = chain_dot(F, self.cov_a, F.T)
return covs
[docs] def cum_errband_mc(self, orth=False, repl=1000,
signif=0.05, seed=None, burn=100):
"""
IRF Monte Carlo integrated error bands of cumulative effect
"""
model = self.model
periods = self.periods
return model.irf_errband_mc(orth=orth, repl=repl,
T=periods, signif=signif, seed=seed, burn=burn, cum=True)
[docs] def lr_effect_cov(self, orth=False):
"""
Returns
-------
"""
lre = self.lr_effects
Finfty = np.kron(np.tile(lre.T, self.lags), lre)
Ik = np.eye(self.neqs)
if orth:
Binf = np.dot(np.kron(self.P.T, np.eye(self.neqs)), Finfty)
Binfbar = np.dot(np.kron(Ik, lre), self.H)
return (chain_dot(Binf, self.cov_a, Binf.T) +
chain_dot(Binfbar, self.cov_sig, Binfbar.T))
else:
return chain_dot(Finfty, self.cov_a, Finfty.T)
[docs] def stderr(self, orth=False):
return np.array([tsa.unvec(np.sqrt(np.diag(c)))
for c in self.cov(orth=orth)])
[docs] def cum_effect_stderr(self, orth=False):
return np.array([tsa.unvec(np.sqrt(np.diag(c)))
for c in self.cum_effect_cov(orth=orth)])
[docs] def lr_effect_stderr(self, orth=False):
cov = self.lr_effect_cov(orth=orth)
return tsa.unvec(np.sqrt(np.diag(cov)))
def _empty_covm(self, periods):
return np.zeros((periods, self.neqs ** 2, self.neqs ** 2),
dtype=float)
[docs] @cache_readonly
def H(self):
k = self.neqs
Lk = tsa.elimination_matrix(k)
Kkk = tsa.commutation_matrix(k, k)
Ik = np.eye(k)
# B = chain_dot(Lk, np.eye(k**2) + commutation_matrix(k, k),
# np.kron(self.P, np.eye(k)), Lk.T)
# return np.dot(Lk.T, L.inv(B))
B = chain_dot(Lk,
np.dot(np.kron(Ik, self.P), Kkk) + np.kron(self.P, Ik),
Lk.T)
return np.dot(Lk.T, L.inv(B))
[docs] def fevd_table(self):
pass
