"""
Vector Autoregression (VAR) processes
References
----------
Lütkepohl (2005) New Introduction to Multiple Time Series Analysis
"""
from statsmodels.compat.pandas import deprecate_kwarg
import numpy as np
import numpy.linalg as npl
from numpy.linalg import slogdet
from statsmodels.tools.numdiff import approx_fprime, approx_hess
from statsmodels.tools.rng_qrng import check_random_state
from statsmodels.tools.validation import string_like
import statsmodels.tsa.base.tsa_model as tsbase
from statsmodels.tsa.vector_ar import util
from statsmodels.tsa.vector_ar.hypothesis_test_results import ErrorBand
from statsmodels.tsa.vector_ar.irf import IRAnalysis
from statsmodels.tsa.vector_ar.var_model import VARProcess, VARResults
def svar_ckerr(svar_type, A, B):
if A is None and (svar_type == "A" or svar_type == "AB"):
raise ValueError("SVAR of type A or AB but A array not given.")
if B is None and (svar_type == "B" or svar_type == "AB"):
raise ValueError("SVAR of type B or AB but B array not given.")
[docs]
class SVAR(tsbase.TimeSeriesModel):
r"""
Fit VAR and then estimate structural components of A and B, defined:
.. math:: Ay_t = A_1 y_{t-1} + \ldots + A_p y_{t-p} + B \varepsilon_t
Parameters
----------
endog : array_like
2-d endogenous response variable. The independent variable.
svar_type : {"A", "B", "AB"}
"A" - estimate structural parameters of A matrix, B assumed = I
"B" - estimate structural parameters of B matrix, A assumed = I
"AB" - estimate structural parameters indicated in both A and B matrix
dates : array_like, optional
must match number of rows of endog
freq : str, optional
The frequency of the time-series. A Pandas offset or 'B', 'D', 'W',
'M', 'A', or 'Q'. This is optional if dates are given.
A : ndarray, optional
neqs x neqs with unknown parameters marked with 'E' for estimate
B : ndarray, optional
neqs x neqs with unknown parameters marked with 'E' for estimate
missing : {"none", "drop", "raise"}, optional
Available options are 'none', 'drop', and 'raise'. If 'none', no nan
checking is done. If 'drop', any observations with nans are dropped.
If 'raise', an error is raised. Default is 'none'.
References
----------
Hamilton (1994) Time Series Analysis
"""
def __init__(
self, endog, svar_type, dates=None, freq=None, A=None, B=None, missing="none"
):
super().__init__(endog, None, dates, freq, missing=missing)
# (self.endog, self.names,
# self.dates) = data_util.interpret_data(endog, names, dates)
self.neqs = self.endog.shape[1]
svar_type = string_like(
svar_type, "svar_type", options=("A", "B", "AB"), lower=False
)
self.svar_type = svar_type
svar_ckerr(svar_type, A, B)
self.A_original = A
self.B_original = B
# Set by fit()/_estimate_svar(), once the lag order is known
self.nobs = None
self.sigma_u = None
# initialize A, B as I if not given
# Initialize SVAR masks
if A is None:
A = np.identity(self.neqs)
self.A_mask = A_mask = np.zeros(A.shape, dtype=bool)
else:
A = A.astype("U")
A_mask = np.logical_or(A == "E", A == "e")
self.A_mask = A_mask
if B is None:
B = np.identity(self.neqs)
self.B_mask = B_mask = np.zeros(B.shape, dtype=bool)
else:
B = B.astype("U")
B_mask = np.logical_or(B == "E", B == "e")
self.B_mask = B_mask
# convert A and B to numeric
# TODO: change this when masked support is better or with formula
# integration
Anum = np.zeros(A.shape, dtype=float)
Anum[~A_mask] = A[~A_mask]
Anum[A_mask] = np.nan
self.A = Anum
Bnum = np.zeros(B.shape, dtype=float)
Bnum[~B_mask] = B[~B_mask]
Bnum[B_mask] = np.nan
self.B = Bnum
# LikelihoodModel.__init__(self, endog)
# super().__init__(endog)
[docs]
def fit(
self,
A_guess=None,
B_guess=None,
maxlags=None,
method="ols",
ic=None,
trend="c",
verbose=False,
s_method="mle",
solver="bfgs",
override=False,
maxiter=500,
maxfun=500,
):
"""
Fit the SVAR model and solve for structural parameters
Parameters
----------
A_guess : array_like, optional
A vector of starting values for all parameters to be estimated
in A.
B_guess : array_like, optional
A vector of starting values for all parameters to be estimated
in B.
maxlags : int, optional
Maximum number of lags to check for order selection, defaults to
12 * (nobs/100.)**(1./4), see select_order function
method : {'ols'}, optional
Estimation method to use
ic : {'aic', 'fpe', 'hqic', 'bic', None}, optional
Information criterion to use for VAR order selection.
aic : Akaike
fpe : Final prediction error
hqic : Hannan-Quinn
bic : Bayesian a.k.a. Schwarz
verbose : bool, optional
Print order selection output to the screen
trend : {"c", "ct", "ctt", "n"}, optional
"c" - add constant
"ct" - constant and trend
"ctt" - constant, linear and quadratic trend
"n" - no constant, no trend
Note that these are prepended to the columns of the dataset.
s_method : {'mle'}, optional
Estimation method for structural parameters
solver : {'nm', 'newton', 'bfgs', 'cg', 'ncg', 'powell'}, optional
Solution method
See statsmodels.base for details
override : bool, optional
If True, returns estimates of A and B without checking
order or rank condition
maxiter : int, optional
Number of iterations to perform in solution method
maxfun : int, optional
Number of function evaluations to perform
Returns
-------
est : SVARResults
Notes
-----
Lütkepohl pp. 146-153
Hamilton pp. 324-336
"""
lags = maxlags
if ic is not None:
selections = self.select_order(maxlags=maxlags, verbose=verbose)
if ic not in selections:
raise ValueError(
f"{ic} not recognized, must be among {sorted(selections)}"
)
lags = selections[ic]
if verbose:
print(f"Using {lags:d} based on {ic} criterion")
elif lags is None:
lags = 1
self.nobs = len(self.endog) - lags
# initialize starting parameters
start_params = self._get_init_params(A_guess, B_guess)
return self._estimate_svar(
start_params,
lags,
trend=trend,
solver=solver,
override=override,
maxiter=maxiter,
maxfun=maxfun,
)
def _get_init_params(self, A_guess, B_guess):
"""Returns either the given starting values or .1 if none are given"""
var_type = self.svar_type.lower()
n_masked_a = self.A_mask.sum()
if var_type in ["ab", "a"]:
if A_guess is None:
A_guess = np.array([0.1] * n_masked_a)
elif len(A_guess) != n_masked_a:
msg = "len(A_guess) = %s, there are %s parameters in A"
raise ValueError(msg % (len(A_guess), n_masked_a))
else:
A_guess = []
n_masked_b = self.B_mask.sum()
if var_type in ["ab", "b"]:
if B_guess is None:
B_guess = np.array([0.1] * n_masked_b)
elif len(B_guess) != n_masked_b:
msg = "len(B_guess) = %s, there are %s parameters in B"
raise ValueError(msg % (len(B_guess), n_masked_b))
else:
B_guess = []
return np.r_[A_guess, B_guess]
def _estimate_svar(
self,
start_params,
lags,
maxiter,
maxfun,
trend="c",
solver="nm",
override=False,
):
"""
Estimate the reduced-form VAR and then solve for structural params
Parameters
----------
start_params : array_like
Starting values for the parameters to be estimated.
lags : int
Number of lags in the VAR.
maxiter : int
Number of iterations to perform in solution method.
maxfun : int
Number of function evaluations to perform.
trend : {"n", "c", "ct", "ctt"}, optional
As per above
solver : str, optional
Solution method
override : bool, optional
If True, returns estimates of A and B without checking
order or rank condition
"""
k_trend = util.get_trendorder(trend)
y = self.endog
z = util.get_var_endog(y, lags, trend=trend, has_constant="raise")
y_sample = y[lags:]
# Lutkepohl p75, about 5x faster than stated formula
var_params = np.linalg.lstsq(z, y_sample, rcond=-1)[0]
resid = y_sample - np.dot(z, var_params)
# Unbiased estimate of covariance matrix $\Sigma_u$ of the white noise
# process $u$
# equivalent definition
# .. math:: \frac{1}{T - Kp - 1} Y^\prime (I_T - Z (Z^\prime Z)^{-1}
# Z^\prime) Y
# Ref: Lutkepohl p.75
# df_resid right now is T - Kp - 1, which is a suggested correction
avobs = len(y_sample)
df_resid = avobs - (self.neqs * lags + k_trend)
sse = np.dot(resid.T, resid)
# TODO: should give users the option to use a dof correction or not
omega = sse / df_resid
self.sigma_u = omega
A, B = self._solve_AB(
start_params, override=override, solver=solver, maxiter=maxiter
)
A_mask = self.A_mask
B_mask = self.B_mask
return SVARResults(
y,
z,
var_params,
omega,
lags,
names=self.endog_names,
trend=trend,
dates=self.data.dates,
model=self,
A=A,
B=B,
A_mask=A_mask,
B_mask=B_mask,
)
[docs]
def loglike(self, params):
"""
Loglikelihood for SVAR model
Parameters
----------
params : array_like
The concatenated unknown parameters of the A and B matrices.
Notes
-----
This method assumes that the autoregressive parameters are
first estimated, then likelihood with structural parameters
is estimated
"""
# TODO: this does not look robust if A or B is None
A = self.A
B = self.B
A_mask = self.A_mask
B_mask = self.B_mask
A_len = len(A[A_mask])
B_len = len(B[B_mask])
if A is not None:
A[A_mask] = params[:A_len]
if B is not None:
B[B_mask] = params[A_len : A_len + B_len]
nobs = self.nobs
neqs = self.neqs
sigma_u = self.sigma_u
W = np.dot(npl.inv(B), A)
trc_in = np.dot(np.dot(W.T, W), sigma_u)
sign, b_logdet = slogdet(B**2) # numpy 1.4 compat
b_slogdet = sign * b_logdet
likl = (
-nobs
/ 2.0
* (
neqs * np.log(2 * np.pi)
- np.log(npl.det(A) ** 2)
+ b_slogdet
+ np.trace(trc_in)
)
)
return likl
[docs]
def score(self, AB_mask):
"""
Return the gradient of the loglike at AB_mask
Parameters
----------
AB_mask : ndarray
Concatenated unknown values of the A and B matrices.
Notes
-----
Return numerical gradient
"""
loglike = self.loglike
if AB_mask.ndim > 1:
AB_mask = AB_mask.ravel()
grad = approx_fprime(AB_mask, loglike, epsilon=1e-8)
# workaround shape of grad if only one parameter #9302
if AB_mask.size == 1 and grad.ndim == 2:
grad = grad.ravel()
return grad
[docs]
def hessian(self, AB_mask):
"""
Returns numerical hessian
Parameters
----------
AB_mask : ndarray
Concatenated unknown values of the A and B matrices.
"""
loglike = self.loglike
if AB_mask.ndim > 1:
AB_mask = AB_mask.ravel()
return approx_hess(AB_mask, loglike)
def _solve_AB(self, start_params, maxiter, override=False, solver="bfgs"):
"""
Solves for MLE estimate of structural parameters
Parameters
----------
start_params : array_like
Starting values for the parameters to be estimated.
maxiter : int
The maximum number of iterations.
override : bool, optional
If True, returns estimates of A and B without checking
order or rank condition
solver : str, optional
Solver to be used. The default is 'bfgs'. Other
choices are 'nm' (Nelder-Mead), 'newton' (Newton-Raphson), 'cg'
conjugate, 'ncg' (non-conjugate gradient), and 'powell'.
Returns
-------
A_solve : ndarray
ML solution for the A matrix.
B_solve : ndarray
ML solution for the B matrix.
"""
# TODO: this could stand a refactor
A_mask = self.A_mask
B_mask = self.B_mask
A = self.A
B = self.B
A_len = len(A[A_mask])
A[A_mask] = start_params[:A_len]
B[B_mask] = start_params[A_len:]
if not override:
J = self._compute_J(A, B)
self.check_order(J)
self.check_rank(J)
else: # TODO: change to a warning?
print("Order/rank conditions have not been checked")
if solver == "bfgs":
kwargs = {"gtol": 1e-5}
else:
kwargs = {}
retvals = (
super()
.fit(
start_params=start_params,
method=solver,
maxiter=maxiter,
disp=False,
**kwargs,
)
.params
)
if retvals.ndim > 1:
retvals = retvals.ravel()
A[A_mask] = retvals[:A_len]
B[B_mask] = retvals[A_len:]
return A, B
def _compute_J(self, A_solve, B_solve):
# first compute appropriate duplication matrix
# taken from Magnus and Neudecker (1980),
# "The Elimination Matrix: Some Lemmas and Applications
# the creation of the D_n matrix follows MN (1980) directly,
# while the rest follows Hamilton (1994)
neqs = self.neqs
sigma_u = self.sigma_u
A_mask = self.A_mask
B_mask = self.B_mask
# first generate duplication matrix, see MN (1980) for notation
D_nT = np.zeros([int((1.0 / 2) * (neqs) * (neqs + 1)), neqs**2])
for j in range(neqs):
i = j
while j <= i < neqs:
u = np.zeros([int((1.0 / 2) * neqs * (neqs + 1)), 1])
u[int(j * neqs + (i + 1) - (1.0 / 2) * (j + 1) * j - 1)] = 1
Tij = np.zeros([neqs, neqs])
Tij[i, j] = 1
Tij[j, i] = 1
D_nT = D_nT + np.dot(u, (Tij.ravel("F")[:, None]).T)
i = i + 1
D_n = D_nT.T
D_pl = npl.pinv(D_n)
# generate S_B
S_B = np.zeros((neqs**2, len(A_solve[A_mask])))
S_D = np.zeros((neqs**2, len(B_solve[B_mask])))
j = 0
j_d = 0
if len(A_solve[A_mask]) != 0:
A_vec = np.ravel(A_mask, order="F")
for k in range(neqs**2):
if A_vec[k]:
S_B[k, j] = -1
j += 1
if len(B_solve[B_mask]) != 0:
B_vec = np.ravel(B_mask, order="F")
for k in range(neqs**2):
if B_vec[k]:
S_D[k, j_d] = 1
j_d += 1
# now compute J
invA = npl.inv(A_solve)
J_p1i = np.dot(np.dot(D_pl, np.kron(sigma_u, invA)), S_B)
J_p1 = -2.0 * J_p1i
J_p2 = np.dot(np.dot(D_pl, np.kron(invA, invA)), S_D)
J = np.append(J_p1, J_p2, axis=1)
return J
[docs]
def check_order(self, J):
if np.size(J, axis=0) < np.size(J, axis=1):
raise ValueError("Order condition not met: solution may not be unique")
[docs]
def check_rank(self, J):
rank = np.linalg.matrix_rank(J)
if rank < np.size(J, axis=1):
raise ValueError("Rank condition not met: solution may not be unique.")
[docs]
class SVARProcess(VARProcess):
"""
Class represents a known SVAR(p) process
Parameters
----------
coefs : ndarray (p x k x k)
intercept : ndarray (length k)
sigma_u : ndarray (k x k)
A_solve : ndarray
neqs x neqs array of estimated structural parameters for the A
matrix.
B_solve : ndarray
neqs x neqs array of estimated structural parameters for the B
matrix.
names : sequence of str, optional
length k
"""
def __init__(self, coefs, intercept, sigma_u, A_solve, B_solve, names=None):
self.k_ar = len(coefs)
self.neqs = coefs.shape[1]
self.coefs = coefs
self.intercept = intercept
self.sigma_u = sigma_u
self.A_solve = A_solve
self.B_solve = B_solve
self.names = names
[docs]
def orth_ma_rep(self, maxn=10, P=None):
"""Unavailable for SVAR"""
raise NotImplementedError
[docs]
def svar_ma_rep(self, maxn=10, P=None):
"""
Compute Structural MA coefficient matrices using MLE of A, B
Parameters
----------
maxn : int, optional
Number of MA matrices to compute.
P : ndarray, optional
Solved A^{-1}B coefficient matrix used to orthogonalize the MA
representation. If not provided, it is computed from the
estimated A and B matrices.
Returns
-------
ndarray
The structural MA coefficient matrices.
"""
if P is None:
A_solve = self.A_solve
B_solve = self.B_solve
P = np.dot(npl.inv(A_solve), B_solve)
ma_mats = self.ma_rep(maxn=maxn)
return np.array([np.dot(coefs, P) for coefs in ma_mats])
[docs]
class SVARResults(SVARProcess, VARResults):
"""
Estimate VAR(p) process with fixed number of lags
Parameters
----------
endog : ndarray
endog_lagged : ndarray
params : ndarray
sigma_u : ndarray
lag_order : int
A : ndarray, optional
neqs x neqs array of estimated structural parameters for the A
matrix.
B : ndarray, optional
neqs x neqs array of estimated structural parameters for the B
matrix.
A_mask : ndarray of bool, optional
neqs x neqs array, True where the corresponding entry of `A` is an
estimated (unknown) parameter.
B_mask : ndarray of bool, optional
neqs x neqs array, True where the corresponding entry of `B` is an
estimated (unknown) parameter.
model : SVAR, optional
The SVAR model instance used to produce the results.
trend : {"n", "c", "ct", "ctt"}, optional
names : sequence of str, optional
List of names of the endogenous variables in order of appearance in `endog`.
dates : array_like, optional
Attributes
----------
aic : float
bic : float
bse : ndarray
coefs : ndarray (p x K x K)
Estimated A_i matrices, A_i = coefs[i-1]
cov_params
dates : array_like
detomega : float
df_model : int
df_resid : int
endog : ndarray
endog_lagged : ndarray
fittedvalues : ndarray
fpe : float
intercept : ndarray
info_criteria : dict of str to float
k_ar : int
Order of VAR process
k_trend : int
llf : float
model : SVAR
names : list of str
variable names
neqs : int
Number of variables (equations)
nobs : int
n_totobs : int
params : ndarray ((Kp + 1) x K)
A_i matrices and intercept in stacked form [int A_1 ... A_p]
pvalues : ndarray
resid : ndarray
sigma_u : ndarray (K x K)
Estimate of white noise process variance Var[u_t]
sigma_u_mle : ndarray
stderr : ndarray
trendorder : int or None
tvalues : ndarray
"""
_model_type = "SVAR"
def __init__(
self,
endog,
endog_lagged,
params,
sigma_u,
lag_order,
A=None,
B=None,
A_mask=None,
B_mask=None,
model=None,
trend="c",
names=None,
dates=None,
):
self.model = model
self.endog = endog
self.endog_lagged = endog_lagged
self.dates = dates
self.n_totobs, self.neqs = self.endog.shape
self.nobs = self.n_totobs - lag_order
k_trend = util.get_trendorder(trend)
if k_trend > 0: # make this the polynomial trend order
trendorder = k_trend - 1
else:
trendorder = None
self.k_trend = k_trend
self.k_exog = k_trend # now (0.9) required by VARProcess
self.k_exog_user = 0
self.trendorder = trendorder
self.exog_names = util.make_lag_names(names, lag_order, k_trend)
self.params = params
self.sigma_u = sigma_u
# Each matrix needs to be transposed
reshaped = self.params[self.k_trend :]
reshaped = reshaped.reshape((lag_order, self.neqs, self.neqs))
# Need to transpose each coefficient matrix
intercept = self.params[0]
coefs = reshaped.swapaxes(1, 2).copy()
# SVAR components
# TODO: if you define these here, you do not also have to define
# them in SVAR process, but I left them for now -ss
self.A = A
self.B = B
self.A_mask = A_mask
self.B_mask = B_mask
super().__init__(coefs, intercept, sigma_u, A, B, names=names)
[docs]
def irf(self, periods=10, var_order=None):
"""
Analyze structural impulse responses to shocks in system
Parameters
----------
periods : int, optional
var_order : sequence, optional
Alternate variable order for the structural impulse matrix. Not
implemented, and so must be None.
Returns
-------
irf : IRAnalysis
Raises
------
NotImplementedError
If `var_order` is not None.
"""
if var_order is not None:
raise NotImplementedError(
"alternate variable order not implemented (yet)"
)
A = self.A
B = self.B
P = np.dot(npl.inv(A), B)
return IRAnalysis(self, P=P, periods=periods, svar=True)
[docs]
@deprecate_kwarg("seed", "rng")
def sirf_errband_mc(
self,
orth=False,
repl=1000,
steps=10,
signif=0.05,
rng=None,
burn=100,
cum=False,
):
"""
Compute Monte Carlo integrated error bands assuming normally
distributed for impulse response functions
Parameters
----------
orth : bool, optional
Compute orthogonalized impulse response error bands
repl : int, optional
number of Monte Carlo replications to perform
steps : int, optional
number of impulse response periods
signif : float, optional
Significance level for error bars (0 < signif < 1), defaults to
95% CI
rng : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
Source of random numbers used for the Monte Carlo
replications. If `rng` is None, a new ``Generator`` is
created using fresh entropy from the operating system. If
`rng` is an int, a new ``RandomState`` instance is created,
seeded with `rng`; this integer-seeding behavior is
deprecated and will change to creating a ``Generator`` in a
future release. If `rng` is already a ``Generator`` or
``RandomState`` instance, that instance is used.
seed : int, array_like of int, numpy.random.Generator, or numpy.random.RandomState, optional
.. deprecated:: 0.15
seed has been deprecated. In-line with SPEC-007, use
rng for passing a random number generator or seed.
burn : int, optional
number of initial observations to discard for simulation
cum : bool, optional
produce cumulative irf error bands
Returns
-------
ErrorBand
A result object with fields ``lower`` and ``upper``, arrays of
ma_rep Monte Carlo standard errors.
Notes
-----
Lütkepohl (2005) Appendix D
"""
neqs = self.neqs
k_ar = self.k_ar
coefs = self.coefs
sigma_u = self.sigma_u
intercept = self.intercept
nobs = self.nobs
ma_coll = np.zeros((repl, steps + 1, neqs, neqs))
A = self.A
B = self.B
A_mask = self.A_mask
B_mask = self.B_mask
A_pass = self.model.A_original
B_pass = self.model.B_original
s_type = self.model.svar_type
g_list = []
def agg(impulses):
if cum:
return impulses.cumsum(axis=0)
return impulses
# Normalize once so that the same generator instance is advanced
# across replications instead of being re-seeded from the same
# raw seed on every iteration, which would otherwise make every
# replication identical.
rng = check_random_state(rng, deprecated=True)
opt_A = A[A_mask]
opt_B = B[B_mask]
for i in range(repl):
# discard first hundred to correct for starting bias
sim = util.varsim(coefs, intercept, sigma_u, rng=rng, steps=nobs + burn)
sim = sim[burn:]
smod = SVAR(sim, svar_type=s_type, A=A_pass, B=B_pass)
if i == 10:
# Use first 10 to update starting val for remainder of fits
mean_AB = np.mean(g_list, axis=0)
split = len(A[A_mask])
opt_A = mean_AB[:split]
opt_B = mean_AB[split:]
sres = smod.fit(maxlags=k_ar, A_guess=opt_A, B_guess=opt_B)
if i < 10:
# save estimates for starting val if in first 10
g_list.append(
np.append(sres.A[A_mask].tolist(), sres.B[B_mask].tolist())
)
ma_coll[i] = agg(sres.svar_ma_rep(maxn=steps))
ma_sort = np.sort(ma_coll, axis=0) # sort to get quantiles
index = (
int(round(signif / 2 * repl) - 1),
int(round((1 - signif / 2) * repl) - 1),
)
lower = ma_sort[index[0], :, :, :]
upper = ma_sort[index[1], :, :, :]
return ErrorBand(lower, upper)
