Source code for snowdrop.src.numeric.solver.Villemot

# -*- coding: utf-8 -*-
"""
Created on Fri Nov 1, 2019

@author: agoumilevski
"""
from __future__ import division

import numpy as np
import scipy.linalg as la
from sys import exit
from snowdrop.src.misc.termcolor import cprint
from snowdrop.src.numeric.solver.util import getParameters
from snowdrop.src.numeric.solver.util import sorter
from snowdrop.src.preprocessor.function import get_function_and_jacobian
from snowdrop.src.utils.equations import getLeadLagIncidence

count = 1

[docs] def solve_villemot(model,steady_state,p=None,suppress_warnings=False): """ Finds first-order accuracy approximation solution. It implements an algoritm described by Sebastien Villemot in his paper: "Solving rational expectations model at first order: what Dynare does." For references please see https://www.dynare.org/wp-repo/dynarewp002.pdf Parameters: :param model: The Model object. :type model: Instance of a class Model. :param steady_state: Steady state. :type steady_state: list. :param p: Model parameters. :type p: list. :param suppress_warnings: Do not show warnings if True :type suppress_warnings: bool. """ global count if count == 1: count += 1 cprint("Sebastien Villemot solver","blue") print() A,C,R,Z,Xa,J,Ru,U = None,None,None,None,None,None,None,None p = getParameters(parameters=p,model=model) try: #Find jacobian z = np.vstack((steady_state,steady_state,steady_state)) K,jacob = get_function_and_jacobian(model,params=p,y=z,order=1) n = len(jacob) if model.max_lead == 0 and model.min_lag < 0: C = jacob[:,n:2*n] L = jacob[:,2*n:3*n] R = jacob[:,3*n:] C_inv = la.inv(C) A = -C_inv @ L R = -C_inv @ R C = -C_inv @ K model.linear_model["A"] = A model.linear_model["C"] = C model.linear_model["R"] = R model.linear_model["Xa"] = 0*C model.linear_model["Ru"] = 0*R model.linear_model["J"] = 0*A model.linear_model["U"] = 0*A model.linear_model["Z"] = -C_inv else: jac = jacob[:,:3*n] r = jacob[:,3*n:] U,A,C,R,Z,Xa,J,Ru = solution(model=model,Jacob=jac,Psi=-r,K=K,suppress_warnings=suppress_warnings) model.linear_model["A"] = A model.linear_model["C"] = C model.linear_model["R"] = R model.linear_model["Ru"] = Ru model.linear_model["Z"] = np.real(Z) model.linear_model["Xa"] = Xa model.linear_model["J"] = J model.linear_model["U"] = U except : if not suppress_warnings: exit('Villemot solver: unable to find solution of a linear model.\n Please use other solvers: BinderPesaran, AndersonMoore, ...') return model
[docs] def solution(model,Jacob,Psi,K=None,qz_factor=1+1.e-6,debug=False,suppress_warnings=False): """ Find first-order accuracy approximation solution. It implements an algoritm described by Sebastien Villemot in the paper: "Solving rational expectations model at first order: what Dynare does." For references please see https://www.dynare.org/wp-repo/dynarewp002.pdf Parameters: :param model: The Model object. :type model: instance of class Model. :param jacob: Jacobian matrix. :type jacob: numpy.ndarray. :param Psi: Matrix of coefficients of shocks. :type Psi: numpy.ndarray. :param C: Array of constants. :type C: numpy.ndarray. """ debug &= suppress_warnings variables = model.symbols["variables"] n = len(Jacob) n_shock = Psi.shape[1] # n_shocks = Psi.shape[1] A,B,R = None,None,None if K is None: K = np.zeros(n) # Retrieve lead-lag incidence matrix lli = model.lead_lag_incidence if lli is None: lli = getLeadLagIncidence(model) forward_columns = np.array([int(f) for f in lli[2] if not np.isnan(f)]) forward_variables = [x for x in model.topology["forward"].keys() if not "_p_" in x] forward_variables = [variables[i] for i in forward_columns] N_forward = len(forward_columns) N = n + N_forward ### ------------------- Make structural state-space repesentation of equations Psi1 = np.concatenate((Psi,np.zeros((N_forward,n_shock)))) K1 = np.concatenate((K,np.zeros(N_forward))) # Build matrix A A = np.zeros((N,N)) # Current values A[:n,:n] = Jacob[:,n:2*n] # Lead variables values A[:n,n:] = Jacob[:,forward_columns] # Identity equations for i in range(N_forward): col = forward_columns[i] A[n+i,col] = 1 # Build matrix B B = np.zeros((N,N)) # Lag variables values B[:n,:n] = -Jacob[:,2*n:3*n] # Identity equations B[n:,n:] = np.eye(N_forward) if debug: np.set_printoptions(precision=1) print("Matrix A:") print(A) print() print("Matrix B:") print(B) print() ### ------------------------ Generalized Schur decomposition # Apply generalized Schur decomposition (a.k.a. QZ decomposition) to the pencil (A,B) a,b,alpha,beta,q,z = la.ordqz(A,B,output='complex',sort=sorter) # Transpose matrices since python QZ algorithm returns a transposed matrix versus Matlab z = z.T.conj() # Check correctness of QZ decomposition err1 = la.norm(q @ a @ z - A) / la.norm(A) err2 = la.norm(q @ b @ z - B) / la.norm(B) if err1 > 1.e-10 or err2 > 1.e-10: cprint("Villemot Solver: QZ decomposition error. \n Inconsistency of A and B matrix decomposition of {0} and {1}".format(round(err1,4),round(err2,4)),"red") raise # Find number of stable and unstable roots roots = np.array([abs(y/x) if abs(x)>0 else np.inf for x,y in zip(np.diag(a),np.diag(b)) ]) n_unstable = sum(roots>=qz_factor) n_stable = N - n_unstable # Check Blanchard-Kahn condition if not suppress_warnings: print() if N_forward == n_unstable: cprint('Blanchard-Kahn condition of existance and unique solution is satisfied','green') elif N_forward < n_unstable: raise Exception('Blanchard-Kahn condition is not satisfied: no stable solution!') else: raise Exception('Blanchard-Kahn condition is not satisfied: multiple stable solutions!') print() if debug: print('n_unstable=',n_unstable,', n_forward_looking=',N_forward,' forward variables=',forward_variables) print('roots:',roots) unstable_roots = roots[n_stable:] if np.any(unstable_roots<1) and not suppress_warnings: cprint("Some of unstable roots are less than one: \n {}".format(np.round(unstable_roots,2)),"red") ### --------------------------Find solution for dynamic variables ### Partition matrices a11 = a[:n_stable, :n_stable] a12 = a[:n_stable, n_stable:] a22 = a[n_stable:, n_stable:] b11 = b[:n_stable, :n_stable] b12 = b[:n_stable, n_stable:] b22 = b[n_stable:, n_stable:] z11 = z[:n_stable,:n_stable] z12 = z[:n_stable,n_stable:] z21 = z[n_stable:,:n_stable] z22 = z[n_stable:,n_stable:] U = z11 D = q @ Psi1 #psi1 = D[:n_stable] psi2 = D[n_stable:] #Const = q @ K1 #k1 = Const[:n_stable] #k2 = Const[n_stable:] # Non-predetermined (or unstable) variables solution try: g_plus = -la.solve(z22, z21) except la.LinAlgError as err: if not suppress_warnings: cprint(err,"red") cprint("n={0}, n_unstable={1}".format(n,n_unstable),"red") print('det(z22)=',la.det(z22)) raise # Predetermined (state) variables solution X = z11 + z12 @ g_plus try: tmp = la.solve(a11,b11 @ X) except la.LinAlgError as err: if not suppress_warnings: cprint(err,"red") cprint("n={0}, n_unstable={1}".format(n,n_unstable),"red") print('det(b11)=',la.det(b11)) raise try: g_minus = la.solve(X,tmp) except la.LinAlgError as err: if not suppress_warnings: cprint(err,"red") cprint("n={0}, n_unstable={1}".format(n,n_unstable),"red") print('det(X)=',la.det(X)) raise ### ----------------------------- Build matrices # Build transition matrix A = np.real(g_minus) # Build matrix R of shocks Fy0 = Jacob[:,n:2*n] Fyp = Jacob[:,forward_columns] temp = Fyp @ g_plus + Fy0 try: Z = la.inv(temp) except la.LinAlgError as err: if not suppress_warnings: cprint(err.message,"red") print('Z=',temp) raise # Shock matrix R = Z @ Psi # Constant matrix C = -np.real(Z @ K) # Get solution decomposition matrices for future shocks # Unstable block G = -la.solve(z11,z12) Ru = -la.solve(b22,psi2) Xa0 = la.solve(a11,b11 @ G + b12) Xa1 = G + la.solve(a11,a12) J = -la.solve(b22,a22) Xa = Xa1 + Xa0 @ J #Xa = G @ J - la.inv(a11) @ b12 return U,A,C,R,Z,Xa,J,Ru