1# -*- coding: utf-8 -*- 2""" 3Impulse reponse-related code 4""" 5 6import numpy as np 7import numpy.linalg as la 8import scipy.linalg as L 9 10from statsmodels.tools.decorators import cache_readonly 11import statsmodels.tsa.tsatools as tsa 12import statsmodels.tsa.vector_ar.plotting as plotting 13import statsmodels.tsa.vector_ar.util as util 14 15mat = np.array 16 17 18class BaseIRAnalysis(object): 19 """ 20 Base class for plotting and computing IRF-related statistics, want to be 21 able to handle known and estimated processes 22 """ 23 24 def __init__(self, model, P=None, periods=10, order=None, svar=False, 25 vecm=False): 26 self.model = model 27 self.periods = periods 28 self.neqs, self.lags, self.T = model.neqs, model.k_ar, model.nobs 29 30 self.order = order 31 32 if P is None: 33 sigma = model.sigma_u 34 35 # TODO, may be difficult at the moment 36 # if order is not None: 37 # indexer = [model.get_eq_index(name) for name in order] 38 # sigma = sigma[:, indexer][indexer, :] 39 40 # if sigma.shape != model.sigma_u.shape: 41 # raise ValueError('variable order is wrong length') 42 43 P = la.cholesky(sigma) 44 45 self.P = P 46 47 self.svar = svar 48 49 self.irfs = model.ma_rep(periods) 50 if svar: 51 self.svar_irfs = model.svar_ma_rep(periods, P=P) 52 else: 53 self.orth_irfs = model.orth_ma_rep(periods, P=P) 54 55 self.cum_effects = self.irfs.cumsum(axis=0) 56 if svar: 57 self.svar_cum_effects = self.svar_irfs.cumsum(axis=0) 58 else: 59 self.orth_cum_effects = self.orth_irfs.cumsum(axis=0) 60 61 # long-run effects may be infinite for VECMs. 62 if not vecm: 63 self.lr_effects = model.long_run_effects() 64 if svar: 65 self.svar_lr_effects = np.dot(model.long_run_effects(), P) 66 else: 67 self.orth_lr_effects = np.dot(model.long_run_effects(), P) 68 69 # auxiliary stuff 70 if vecm: 71 self._A = util.comp_matrix(model.var_rep) 72 else: 73 self._A = util.comp_matrix(model.coefs) 74 75 def _choose_irfs(self, orth=False, svar=False): 76 if orth: 77 return self.orth_irfs 78 elif svar: 79 return self.svar_irfs 80 else: 81 return self.irfs 82 83 def cov(self, *args, **kwargs): 84 raise NotImplementedError 85 86 def cum_effect_cov(self, *args, **kwargs): 87 raise NotImplementedError 88 89 def plot(self, orth=False, *, impulse=None, response=None, 90 signif=0.05, plot_params=None, figsize=(10, 10), 91 subplot_params=None, plot_stderr=True, stderr_type='asym', 92 repl=1000, seed=None, component=None): 93 """ 94 Plot impulse responses 95 96 Parameters 97 ---------- 98 orth : bool, default False 99 Compute orthogonalized impulse responses 100 impulse : {str, int} 101 variable providing the impulse 102 response : {str, int} 103 variable affected by the impulse 104 signif : float (0 < signif < 1) 105 Significance level for error bars, defaults to 95% CI 106 subplot_params : dict 107 To pass to subplot plotting funcions. Example: if fonts are too big, 108 pass {'fontsize' : 8} or some number to your taste. 109 plot_params : dict 110 111 figsize : (float, float), default (10, 10) 112 Figure size (width, height in inches) 113 plot_stderr : bool, default True 114 Plot standard impulse response error bands 115 stderr_type : str 116 'asym': default, computes asymptotic standard errors 117 'mc': monte carlo standard errors (use rpl) 118 repl : int, default 1000 119 Number of replications for Monte Carlo and Sims-Zha standard errors 120 seed : int 121 np.random.seed for Monte Carlo replications 122 component: array or vector of principal component indices 123 """ 124 periods = self.periods 125 model = self.model 126 svar = self.svar 127 128 if orth and svar: 129 raise ValueError("For SVAR system, set orth=False") 130 131 irfs = self._choose_irfs(orth, svar) 132 if orth: 133 title = 'Impulse responses (orthogonalized)' 134 elif svar: 135 title = 'Impulse responses (structural)' 136 else: 137 title = 'Impulse responses' 138 139 if plot_stderr is False: 140 stderr = None 141 142 elif stderr_type not in ['asym', 'mc', 'sz1', 'sz2','sz3']: 143 raise ValueError("Error type must be either 'asym', 'mc','sz1','sz2', or 'sz3'") 144 else: 145 if stderr_type == 'asym': 146 stderr = self.cov(orth=orth) 147 if stderr_type == 'mc': 148 stderr = self.errband_mc(orth=orth, svar=svar, 149 repl=repl, signif=signif, 150 seed=seed) 151 if stderr_type == 'sz1': 152 stderr = self.err_band_sz1(orth=orth, svar=svar, 153 repl=repl, signif=signif, 154 seed=seed, 155 component=component) 156 if stderr_type == 'sz2': 157 stderr = self.err_band_sz2(orth=orth, svar=svar, 158 repl=repl, signif=signif, 159 seed=seed, 160 component=component) 161 if stderr_type == 'sz3': 162 stderr = self.err_band_sz3(orth=orth, svar=svar, 163 repl=repl, signif=signif, 164 seed=seed, 165 component=component) 166 167 fig = plotting.irf_grid_plot(irfs, stderr, impulse, response, 168 self.model.names, title, signif=signif, 169 subplot_params=subplot_params, 170 plot_params=plot_params, 171 figsize=figsize, 172 stderr_type=stderr_type) 173 return fig 174 175 def plot_cum_effects(self, orth=False, *, impulse=None, response=None, 176 signif=0.05, plot_params=None, figsize=(10, 10), 177 subplot_params=None, plot_stderr=True, 178 stderr_type='asym', repl=1000, seed=None): 179 """ 180 Plot cumulative impulse response functions 181 182 Parameters 183 ---------- 184 orth : bool, default False 185 Compute orthogonalized impulse responses 186 impulse : {str, int} 187 variable providing the impulse 188 response : {str, int} 189 variable affected by the impulse 190 signif : float (0 < signif < 1) 191 Significance level for error bars, defaults to 95% CI 192 subplot_params : dict 193 To pass to subplot plotting funcions. Example: if fonts are too big, 194 pass {'fontsize' : 8} or some number to your taste. 195 plot_params : dict 196 197 figsize: (float, float), default (10, 10) 198 Figure size (width, height in inches) 199 plot_stderr : bool, default True 200 Plot standard impulse response error bands 201 stderr_type : str 202 'asym': default, computes asymptotic standard errors 203 'mc': monte carlo standard errors (use rpl) 204 repl : int, default 1000 205 Number of replications for monte carlo standard errors 206 seed : int 207 np.random.seed for Monte Carlo replications 208 """ 209 210 if orth: 211 title = 'Cumulative responses responses (orthogonalized)' 212 cum_effects = self.orth_cum_effects 213 lr_effects = self.orth_lr_effects 214 else: 215 title = 'Cumulative responses' 216 cum_effects = self.cum_effects 217 lr_effects = self.lr_effects 218 219 if stderr_type not in ['asym', 'mc']: 220 raise ValueError("`stderr_type` must be one of 'asym', 'mc'") 221 else: 222 if stderr_type == 'asym': 223 stderr = self.cum_effect_cov(orth=orth) 224 if stderr_type == 'mc': 225 stderr = self.cum_errband_mc(orth=orth, repl=repl, 226 signif=signif, seed=seed) 227 if not plot_stderr: 228 stderr = None 229 230 fig = plotting.irf_grid_plot(cum_effects, stderr, impulse, response, 231 self.model.names, title, signif=signif, 232 hlines=lr_effects, 233 subplot_params=subplot_params, 234 plot_params=plot_params, 235 figsize=figsize, 236 stderr_type=stderr_type) 237 return fig 238 239 240class IRAnalysis(BaseIRAnalysis): 241 """ 242 Impulse response analysis class. Computes impulse responses, asymptotic 243 standard errors, and produces relevant plots 244 245 Parameters 246 ---------- 247 model : VAR instance 248 249 Notes 250 ----- 251 Using Lütkepohl (2005) notation 252 """ 253 def __init__(self, model, P=None, periods=10, order=None, svar=False, 254 vecm=False): 255 BaseIRAnalysis.__init__(self, model, P=P, periods=periods, 256 order=order, svar=svar, vecm=vecm) 257 258 if vecm: 259 self.cov_a = model.cov_var_repr 260 else: 261 self.cov_a = model._cov_alpha 262 self.cov_sig = model._cov_sigma 263 264 # memoize dict for G matrix function 265 self._g_memo = {} 266 267 def cov(self, orth=False): 268 """ 269 Compute asymptotic standard errors for impulse response coefficients 270 271 Notes 272 ----- 273 Lütkepohl eq 3.7.5 274 275 Returns 276 ------- 277 """ 278 if orth: 279 return self._orth_cov() 280 281 covs = self._empty_covm(self.periods + 1) 282 covs[0] = np.zeros((self.neqs ** 2, self.neqs ** 2)) 283 for i in range(1, self.periods + 1): 284 Gi = self.G[i - 1] 285 covs[i] = Gi @ self.cov_a @ Gi.T 286 287 return covs 288 289 def errband_mc(self, orth=False, svar=False, repl=1000, 290 signif=0.05, seed=None, burn=100): 291 """ 292 IRF Monte Carlo integrated error bands 293 """ 294 model = self.model 295 periods = self.periods 296 if svar: 297 return model.sirf_errband_mc(orth=orth, repl=repl, steps=periods, 298 signif=signif, seed=seed, 299 burn=burn, cum=False) 300 else: 301 return model.irf_errband_mc(orth=orth, repl=repl, steps=periods, 302 signif=signif, seed=seed, 303 burn=burn, cum=False) 304 305 def err_band_sz1(self, orth=False, svar=False, repl=1000, 306 signif=0.05, seed=None, burn=100, component=None): 307 """ 308 IRF Sims-Zha error band method 1. Assumes symmetric error bands around 309 mean. 310 311 Parameters 312 ---------- 313 orth : bool, default False 314 Compute orthogonalized impulse responses 315 repl : int, default 1000 316 Number of MC replications 317 signif : float (0 < signif < 1) 318 Significance level for error bars, defaults to 95% CI 319 seed : int, default None 320 np.random seed 321 burn : int, default 100 322 Number of initial simulated obs to discard 323 component : neqs x neqs array, default to largest for each 324 Index of column of eigenvector/value to use for each error band 325 Note: period of impulse (t=0) is not included when computing 326 principle component 327 328 References 329 ---------- 330 Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse 331 Response". Econometrica 67: 1113-1155. 332 """ 333 334 model = self.model 335 periods = self.periods 336 irfs = self._choose_irfs(orth, svar) 337 neqs = self.neqs 338 irf_resim = model.irf_resim(orth=orth, repl=repl, steps=periods, 339 seed=seed, burn=burn) 340 q = util.norm_signif_level(signif) 341 342 W, eigva, k =self._eigval_decomp_SZ(irf_resim) 343 344 if component is not None: 345 if np.shape(component) != (neqs,neqs): 346 raise ValueError("Component array must be " + str(neqs) + " x " + str(neqs)) 347 if np.argmax(component) >= neqs*periods: 348 raise ValueError("Atleast one of the components does not exist") 349 else: 350 k = component 351 352 # here take the kth column of W, which we determine by finding the largest eigenvalue of the covaraince matrix 353 lower = np.copy(irfs) 354 upper = np.copy(irfs) 355 for i in range(neqs): 356 for j in range(neqs): 357 lower[1:,i,j] = irfs[1:,i,j] + W[i,j,:,k[i,j]]*q*np.sqrt(eigva[i,j,k[i,j]]) 358 upper[1:,i,j] = irfs[1:,i,j] - W[i,j,:,k[i,j]]*q*np.sqrt(eigva[i,j,k[i,j]]) 359 360 return lower, upper 361 362 def err_band_sz2(self, orth=False, svar=False, repl=1000, signif=0.05, 363 seed=None, burn=100, component=None): 364 """ 365 IRF Sims-Zha error band method 2. 366 367 This method Does not assume symmetric error bands around mean. 368 369 Parameters 370 ---------- 371 orth : bool, default False 372 Compute orthogonalized impulse responses 373 repl : int, default 1000 374 Number of MC replications 375 signif : float (0 < signif < 1) 376 Significance level for error bars, defaults to 95% CI 377 seed : int, default None 378 np.random seed 379 burn : int, default 100 380 Number of initial simulated obs to discard 381 component : neqs x neqs array, default to largest for each 382 Index of column of eigenvector/value to use for each error band 383 Note: period of impulse (t=0) is not included when computing 384 principle component 385 386 References 387 ---------- 388 Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse 389 Response". Econometrica 67: 1113-1155. 390 """ 391 model = self.model 392 periods = self.periods 393 irfs = self._choose_irfs(orth, svar) 394 neqs = self.neqs 395 irf_resim = model.irf_resim(orth=orth, repl=repl, steps=periods, seed=seed, 396 burn=100) 397 398 W, eigva, k = self._eigval_decomp_SZ(irf_resim) 399 400 if component is not None: 401 if np.shape(component) != (neqs,neqs): 402 raise ValueError("Component array must be " + str(neqs) + " x " + str(neqs)) 403 if np.argmax(component) >= neqs*periods: 404 raise ValueError("Atleast one of the components does not exist") 405 else: 406 k = component 407 408 gamma = np.zeros((repl, periods+1, neqs, neqs)) 409 for p in range(repl): 410 for i in range(neqs): 411 for j in range(neqs): 412 gamma[p,1:,i,j] = W[i,j,k[i,j],:] * irf_resim[p,1:,i,j] 413 414 gamma_sort = np.sort(gamma, axis=0) #sort to get quantiles 415 indx = round(signif/2*repl)-1,round((1-signif/2)*repl)-1 416 417 lower = np.copy(irfs) 418 upper = np.copy(irfs) 419 for i in range(neqs): 420 for j in range(neqs): 421 lower[:,i,j] = irfs[:,i,j] + gamma_sort[indx[0],:,i,j] 422 upper[:,i,j] = irfs[:,i,j] + gamma_sort[indx[1],:,i,j] 423 424 return lower, upper 425 426 def err_band_sz3(self, orth=False, svar=False, repl=1000, signif=0.05, 427 seed=None, burn=100, component=None): 428 """ 429 IRF Sims-Zha error band method 3. Does not assume symmetric error bands around mean. 430 431 Parameters 432 ---------- 433 orth : bool, default False 434 Compute orthogonalized impulse responses 435 repl : int, default 1000 436 Number of MC replications 437 signif : float (0 < signif < 1) 438 Significance level for error bars, defaults to 95% CI 439 seed : int, default None 440 np.random seed 441 burn : int, default 100 442 Number of initial simulated obs to discard 443 component : vector length neqs, default to largest for each 444 Index of column of eigenvector/value to use for each error band 445 Note: period of impulse (t=0) is not included when computing 446 principle component 447 448 References 449 ---------- 450 Sims, Christopher A., and Tao Zha. 1999. "Error Bands for Impulse 451 Response". Econometrica 67: 1113-1155. 452 """ 453 454 model = self.model 455 periods = self.periods 456 irfs = self._choose_irfs(orth, svar) 457 neqs = self.neqs 458 irf_resim = model.irf_resim(orth=orth, repl=repl, steps=periods, 459 seed=seed, burn=100) 460 stack = np.zeros((neqs, repl, periods*neqs)) 461 462 #stack left to right, up and down 463 464 for p in range(repl): 465 for i in range(neqs): 466 stack[i, p,:] = np.ravel(irf_resim[p,1:,:,i].T) 467 468 stack_cov=np.zeros((neqs, periods*neqs, periods*neqs)) 469 W = np.zeros((neqs, periods*neqs, periods*neqs)) 470 eigva = np.zeros((neqs, periods*neqs)) 471 k = np.zeros(neqs, dtype=int) 472 473 if component is not None: 474 if np.size(component) != (neqs): 475 raise ValueError("Component array must be of length " + str(neqs)) 476 if np.argmax(component) >= neqs*periods: 477 raise ValueError("Atleast one of the components does not exist") 478 else: 479 k = component 480 481 #compute for eigen decomp for each stack 482 for i in range(neqs): 483 stack_cov[i] = np.cov(stack[i],rowvar=0) 484 W[i], eigva[i], k[i] = util.eigval_decomp(stack_cov[i]) 485 486 gamma = np.zeros((repl, periods+1, neqs, neqs)) 487 for p in range(repl): 488 c = 0 489 for j in range(neqs): 490 for i in range(neqs): 491 gamma[p,1:,i,j] = W[j,k[j],i*periods:(i+1)*periods] * irf_resim[p,1:,i,j] 492 if i == neqs-1: 493 gamma[p,1:,i,j] = W[j,k[j],i*periods:] * irf_resim[p,1:,i,j] 494 495 gamma_sort = np.sort(gamma, axis=0) #sort to get quantiles 496 indx = round(signif/2*repl)-1,round((1-signif/2)*repl)-1 497 498 lower = np.copy(irfs) 499 upper = np.copy(irfs) 500 for i in range(neqs): 501 for j in range(neqs): 502 lower[:,i,j] = irfs[:,i,j] + gamma_sort[indx[0],:,i,j] 503 upper[:,i,j] = irfs[:,i,j] + gamma_sort[indx[1],:,i,j] 504 505 return lower, upper 506 507 def _eigval_decomp_SZ(self, irf_resim): 508 """ 509 Returns 510 ------- 511 W: array of eigenvectors 512 eigva: list of eigenvalues 513 k: matrix indicating column # of largest eigenvalue for each c_i,j 514 """ 515 neqs = self.neqs 516 periods = self.periods 517 518 cov_hold = np.zeros((neqs, neqs, periods, periods)) 519 for i in range(neqs): 520 for j in range(neqs): 521 cov_hold[i,j,:,:] = np.cov(irf_resim[:,1:,i,j],rowvar=0) 522 523 W = np.zeros((neqs, neqs, periods, periods)) 524 eigva = np.zeros((neqs, neqs, periods, 1)) 525 k = np.zeros((neqs, neqs), dtype=int) 526 527 for i in range(neqs): 528 for j in range(neqs): 529 W[i,j,:,:], eigva[i,j,:,0], k[i,j] = util.eigval_decomp(cov_hold[i,j,:,:]) 530 return W, eigva, k 531 532 @cache_readonly 533 def G(self): 534 # Gi matrices as defined on p. 111 535 536 K = self.neqs 537 538 # nlags = self.model.p 539 # J = np.hstack((np.eye(K),) + (np.zeros((K, K)),) * (nlags - 1)) 540 541 def _make_g(i): 542 # p. 111 Lutkepohl 543 G = 0. 544 for m in range(i): 545 # be a bit cute to go faster 546 idx = i - 1 - m 547 if idx in self._g_memo: 548 apow = self._g_memo[idx] 549 else: 550 apow = la.matrix_power(self._A.T, idx) 551 # apow = np.dot(J, apow) 552 apow = apow[:K] 553 self._g_memo[idx] = apow 554 555 # take first K rows 556 piece = np.kron(apow, self.irfs[m]) 557 G = G + piece 558 559 return G 560 561 return [_make_g(i) for i in range(1, self.periods + 1)] 562 563 def _orth_cov(self): 564 # Lutkepohl 3.7.8 565 566 Ik = np.eye(self.neqs) 567 PIk = np.kron(self.P.T, Ik) 568 H = self.H 569 570 covs = self._empty_covm(self.periods + 1) 571 for i in range(self.periods + 1): 572 if i == 0: 573 apiece = 0 574 else: 575 Ci = np.dot(PIk, self.G[i-1]) 576 apiece = Ci @ self.cov_a @ Ci.T 577 578 Cibar = np.dot(np.kron(Ik, self.irfs[i]), H) 579 bpiece = (Cibar @ self.cov_sig @ Cibar.T) / self.T 580 581 # Lutkepohl typo, cov_sig correct 582 covs[i] = apiece + bpiece 583 584 return covs 585 586 def cum_effect_cov(self, orth=False): 587 """ 588 Compute asymptotic standard errors for cumulative impulse response 589 coefficients 590 591 Parameters 592 ---------- 593 orth : bool 594 595 Notes 596 ----- 597 eq. 3.7.7 (non-orth), 3.7.10 (orth) 598 599 Returns 600 ------- 601 """ 602 Ik = np.eye(self.neqs) 603 PIk = np.kron(self.P.T, Ik) 604 605 F = 0. 606 covs = self._empty_covm(self.periods + 1) 607 for i in range(self.periods + 1): 608 if i > 0: 609 F = F + self.G[i - 1] 610 611 if orth: 612 if i == 0: 613 apiece = 0 614 else: 615 Bn = np.dot(PIk, F) 616 apiece = Bn @ self.cov_a @ Bn.T 617 618 Bnbar = np.dot(np.kron(Ik, self.cum_effects[i]), self.H) 619 bpiece = (Bnbar @ self.cov_sig @ Bnbar.T) / self.T 620 621 covs[i] = apiece + bpiece 622 else: 623 if i == 0: 624 covs[i] = np.zeros((self.neqs**2, self.neqs**2)) 625 continue 626 627 covs[i] = F @ self.cov_a @ F.T 628 629 return covs 630 631 def cum_errband_mc(self, orth=False, repl=1000, 632 signif=0.05, seed=None, burn=100): 633 """ 634 IRF Monte Carlo integrated error bands of cumulative effect 635 """ 636 model = self.model 637 periods = self.periods 638 return model.irf_errband_mc(orth=orth, repl=repl, 639 steps=periods, signif=signif, 640 seed=seed, burn=burn, cum=True) 641 642 def lr_effect_cov(self, orth=False): 643 """ 644 Returns 645 ------- 646 """ 647 lre = self.lr_effects 648 Finfty = np.kron(np.tile(lre.T, self.lags), lre) 649 Ik = np.eye(self.neqs) 650 651 if orth: 652 Binf = np.dot(np.kron(self.P.T, np.eye(self.neqs)), Finfty) 653 Binfbar = np.dot(np.kron(Ik, lre), self.H) 654 655 return (Binf @ self.cov_a @ Binf.T + 656 Binfbar @ self.cov_sig @ Binfbar.T) 657 else: 658 return Finfty @ self.cov_a @ Finfty.T 659 660 def stderr(self, orth=False): 661 return np.array([tsa.unvec(np.sqrt(np.diag(c))) 662 for c in self.cov(orth=orth)]) 663 664 def cum_effect_stderr(self, orth=False): 665 return np.array([tsa.unvec(np.sqrt(np.diag(c))) 666 for c in self.cum_effect_cov(orth=orth)]) 667 668 def lr_effect_stderr(self, orth=False): 669 cov = self.lr_effect_cov(orth=orth) 670 return tsa.unvec(np.sqrt(np.diag(cov))) 671 672 def _empty_covm(self, periods): 673 return np.zeros((periods, self.neqs ** 2, self.neqs ** 2), 674 dtype=float) 675 676 @cache_readonly 677 def H(self): 678 k = self.neqs 679 Lk = tsa.elimination_matrix(k) 680 Kkk = tsa.commutation_matrix(k, k) 681 Ik = np.eye(k) 682 683 # B = Lk @ (np.eye(k**2) + commutation_matrix(k, k)) @ \ 684 # np.kron(self.P, np.eye(k)) @ Lk.T 685 # return Lk.T @ L.inv(B) 686 687 B = Lk @ (np.kron(Ik, self.P) @ Kkk + np.kron(self.P, Ik)) @ Lk.T 688 689 return np.dot(Lk.T, L.inv(B)) 690 691 def fevd_table(self): 692 raise NotImplementedError 693