Temporary commit after much work on abcs_rbcs model branch and implementing test cases

Author: escheffel

pymaclab.wpu

@@ -109,7 +109,7 @@ def init1(self,no_wrap=False):
         if txtpars.fmt == 'pymaclab':
             other = populate_model_stage_one(other, secs)
         elif txtpars.fmt == 'dynarepp':
-            pml_modfile = dynarepp_to_pml_translate(secs=secs)
+            pml_modfile = dynarepp_to_pml.translate(secs=secs)
 
         if not no_wrap:
             # Open updaters path

pymaclab/dsge/inits/_inits.py

@@ -68,24 +68,8 @@ def populate_model_stage_one(self, secs):
         audic['iid']['mod'] = []
 
         for x in secs['varvec'][0]:
-            if viiexp.search(x):
-                ma = viiexp.search(x)
-                vari = ma.group('vari').strip()
-                viid = ma.group('viid').strip()
-                varn = ma.group('varn').strip()
-                vtype = ma.group('vtype').strip()
-                mods = ma.group('mod')
-                vardic[vtype]['var'].append([vari,varn,viid])
-                if mods != None:
-                    mods = mods.strip()
-                    if ',' in mods:
-                        vardic[vtype]['mod'].append(mods[1:-1].strip().split(','))
-                    else:
-                        vardic[vtype]['mod'].append([mods[1:-1].strip(),])
-                else:
-                    vardic[vtype]['mod'].append([])
 
-            elif vtexp.search(x):
+            if vtexp.search(x):
                 ma = vtexp.search(x)
                 vari = ma.group('vari').strip()
                 varn = ma.group('varn').strip()
@@ -99,7 +83,7 @@ def populate_model_stage_one(self, secs):
                     else:
                         vardic[vtype]['mod'].append([mods[1:-1].strip(),])
                 else:
-                    vardic[vtype]['mod'].append([])
+                    vardic[vtype]['mod'].append([''])
 
             elif voexp.search(x):
                 ma = voexp.search(x)
@@ -115,7 +99,7 @@ def populate_model_stage_one(self, secs):
                     else:
                         vardic['other']['mod'].append([mods[1:-1].strip(),])
                 else:
-                    vardic['other']['mod'].append([])
+                    vardic['other']['mod'].append([''])
 
         self.nendo = len(vardic['endo']['var'])
         self.nexo = len(vardic['exo']['var'])

pymaclab/dsge/parsers/_dsgeparser.py

@@ -16,6 +16,7 @@
 from pymaclab.dsge.helpers import _helpers as HLP
 from scipy.linalg import qz
 import numpy as np
+import scipy
 import pylab as P
 from matplotlib import pyplot as PLT
 import copy as COP
@@ -294,22 +295,91 @@ def mkmats(self):
         # Change NN matrix by swapping signs
         self.N = self.N*(-np.matlib.eye(self.N.shape[0]))
 
-    def solve(self,sol_method='do_QZ',Xi_select='all'):
+    def solve(self,sol_method='do_QZ2',Xi_select='all'):
         self.output = {}
         self.Xi_select = Xi_select
         self.sol_method = sol_method
         CC = self.C
         try:
-            testo = inv(C)
+            testo = inv(CC)
         except:
             print "Error: Matrix C is not of full rank!"
             return False
         if self.sol_method == 'do_PP':
             self.do_PP(self.Xi_select)
         elif self.sol_method == 'do_QZ':
             self.do_QZ(self.Xi_select,Tol=1.0000e-06)
+        elif self.sol_method == 'do_QZ2':
+            self.do_QZ2(self.Xi_select,Tol=1.0000e-06)
         self.do_QRS()
 
+    def do_QZ2(self,Xi_select='all',Tol=1.0000e-06):
+        # Make uhlig matrices locally available for computations
+        AA = self.A
+        BB = self.B
+        CC = self.C
+        DD = self.D
+        FF = self.F
+        GG = self.G
+        HH = self.H
+        JJ = self.J
+        KK = self.K
+        LL = self.L
+        MM = self.M
+        NN = self.N
+        q_expectational_equ = np.shape(FF)[0]
+        m_states = np.shape(FF)[1]
+        l_equ = np.shape(CC)[0]
+        n_endog = np.shape(CC)[1]
+        k_exog = min(np.shape(NN))
+        
+        # Do some test
+        if matrix_rank(CC) < n_endog:
+            print 'uhlerror: Rank(CC) needs to be at least\n\
+            equal to the number of endogenous variables'
+            return False
+        
+        # Set up and find the solution to the matrix quadratic polynomial
+        I1 = MAT.eye(m_states)
+        Z1 = MAT.zeros((m_states,m_states))       
+        invC = inv(CC);
+        psy = FF-JJ*invC*AA;
+        lambdar = JJ*invC*BB-GG+KK*invC*AA;
+        TT = KK*invC*BB-HH;
+        AA1 = MAT.vstack((MAT.hstack((lambdar,TT)),
+                          MAT.hstack((I1,Z1))
+                          ))
+        AA2 = MAT.vstack((MAT.hstack((psy,Z1)),
+                          MAT.hstack((Z1,I1))
+                          ))
+
+        # Find the generalized eigenvalues and eigenvectors
+        # Scipy returns eigval,eigvec, BUT Matlab returns eigvec,eigval !!!
+        [eigval,eigvecl,eigvecr] = scipy.linalg.eig(AA1,AA2,left=True,right=True)
+        eigval = MAT.real(eigval)
+        eigvecr = MAT.real(eigvecr)
+        
+        # Generate selection index
+        dtype = [('eigvals',float),('index', float)]
+        values = [(x,i1) for i1,x in enumerate(eigval)]
+        eigval = np.array(values,dtype=dtype)
+        eigval = np.sort(eigval,axis=0,order=['eigvals','index'])
+        eigpick = [x for i1,x in enumerate(eigval) if eigval[i1][0] < 1.0]
+        eigindex = np.array([x[1] for x in eigpick])
+        
+        # Use the selection index on eigval
+        eigval = np.array([x[0] for x in eigpick])
+        
+        # Generate DD and EE matrices and calculate PP
+        DD = MAT.zeros([m_states,m_states])
+        for i1 in range(m_states):
+            DD[i1,i1] = eigval[i1]
+        EE = eigvecr[:,list(eigindex)][-m_states:,:]
+        # Find the matrix P
+        PP = EE*DD*inv(EE)
+        self.P = PP
+        self.CC_plus = LIN.pinv(CC)
+        
     def do_QZ(self,Xi_select='all',Tol=1.0000e-06):
         # Make uhlig matrices locally available for computations
         AA = self.A
@@ -371,11 +441,6 @@ def do_QZ(self,Xi_select='all',Tol=1.0000e-06):
                       np.concatenate((MAT.zeros((m_states,m_states)),MAT.identity(m_states)),1)\
                   ))
 
-#        (Delta_up,Xi_up,UUU,VVV)=\
-#         HLP.qz(Delta_mat,Xi_mat,\
-#            mode='complex',\
-#            lapackname=lapackname,\
-#            lapackpath=lapackpath)
         (Delta_up,Xi_up,UUU,VVV) = qz(Delta_mat,Xi_mat)
 
 
@@ -542,16 +607,39 @@ def do_QRS(self):
         (l_equ,m_states) = MAT.shape(AA)
         (l_equ,n_endog) = MAT.shape(CC)
         (l_equ,k_exog) = MAT.shape(DD)
+        
+        # Compute the R matrix
         RR = -CC_plus*(AA*PP+BB)
         self.R = RR
         Qone = MAT.kron(NN.T,(FF-JJ*inv(CC)*AA))
         Qtwo = MAT.kron(MAT.eye(k_exog),(FF*PP+GG+JJ*RR-KK*inv(CC)*AA))
         Q = Qone+Qtwo
         invQ = inv(Q)
-        Qthree = ((JJ*inv(CC)*DD-LL)*NN+KK*inv(CC)*DD-MM).T
-        QQ = (invQ*Qthree).T
-        self.Q = QQ
-        SS = -CC_plus*(AA*QQ+DD)
+        self.invQ = invQ
+        Qthree = ((JJ*inv(CC)*DD-LL)*NN+KK*inv(CC)*DD-MM)
+        self.Qthree = Qthree
+        
+        # Some further manipulation to get Q
+        sizo,sizp = Qthree.shape
+        for i1 in range(sizp):
+            if i1 == 0:
+                Qfour = Qthree[:,i1]
+            else:
+                Qfour = MAT.vstack((Qfour,Qthree[:,i1]))
+        self.Qfour = Qfour
+        QQ = invQ*Qfour
+        self.QQ = QQ
+        for i1 in range(sizp):
+            begini = i1*sizo
+            endi = (i1+1)*sizo
+            if i1 == 0:
+                QQQ = QQ[begini:endi,0]
+            else:
+                QQQ = MAT.hstack((QQQ,QQ[begini:endi,0]))
+        self.Q = QQQ
+        
+        # Compute S matrix
+        SS = -CC_plus*(AA*QQQ+DD)
         self.S = SS
         del self.CC_plus
 

pymaclab/dsge/solvers/modsolvers.py

@@ -1,4 +1,4 @@
-def translate(secs=None,fpath=None):
+def translate(secs=None,fpath=None,dyn_vtimings={'exo':[-1,0],'endo':[-1,0],'iid':[0,1],'con':[0,1]}):
     from pymaclab.modfiles.templates import mako_pml_template
 
     # Render the template to be passed to dynare++

pymaclab/dsge/translators/dynarepp_to_pml.py

@@ -0,0 +1,15 @@
+def translate(template_paramdic=None,fpath=None):
+    import pymaclab.modfiles.templates.wheezy_template as template
+    
+    # Render the template to be passed to wheezy pymaclab template
+    modstr = template.render(template_paramdic)
+    
+    # If a filepath has been passed then just write the pymaclab modfile, but no more!
+    if fpath != None:
+        filo = open(fpath,'w')
+        filo.write(modstr)
+        filo.flush()
+        filo.close()
+        return
+    else:
+        return modstr

pymaclab/dsge/translators/pml_to_pml.py

@@ -1,6 +1,7 @@
 from copy import deepcopy
 from pymaclab.dsge.translators import pml_to_dynarepp
 from pymaclab.dsge.translators import dynarepp_to_pml
+from pymaclab.dsge.translators import pml_to_pml
 
 
 class Translators(object):
@@ -29,4 +30,16 @@ def dynarepp_to_pml(self,template_paramdic=None,fpath=None,focli=None):
             if template_paramdic == None:
                 dynarepp_to_pml.translate(template_paramdic=self.template_paramdic,fpath=fpath,focli=focli)
             else:
-                dynarepp_to_pml.translate(template_paramdic=template_paramdic,fpath=fpath,focli=focli)
+                dynarepp_to_pml.translate(template_paramdic=template_paramdic,fpath=fpath,focli=focli)
+
+    def pml_to_pml(self,template_paramdic=None,fpath=None):
+        if fpath == None:
+            if template_paramdic == None:
+                return pml_to_pml.translate(template_paramdic=self.template_paramdic)
+            else:
+                return pml_to_pml.translate(template_paramdic=template_paramdic)
+        else:
+            if template_paramdic == None:
+                pml_to_pml.translate(template_paramdic=self.template_paramdic,fpath=fpath)
+            else:
+                pml_to_pml.translate(template_paramdic=template_paramdic,fpath=fpath)

pymaclab/dsge/translators/translators.py

@@ -27,15 +27,14 @@ sigma_epsu = 0.01;
 
 %Variable Vectors+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 [1]  k(t):capital{endo}[log,bk]
-[2]  c(t):consumption{con}[log,bk]
-[3]  h(t):labour{con}[log,bk]
-[4]  w(t):real_wage{con}[log,bk]
-[5]  r(t):real_rate_capital{con}[log,bk]
-[6]  p(t):rprice{con}[log,bk]
+[4]  r(t):real_rate_capital{con}[log,bk]
+[3]  w(t):real_wage{con}[log,bk]
+[2]  h(t):labour{con}[log,bk]
+[5]  p(t):rprice{con}[log,bk]
+[6]  lam(t):productivity{exo}[log,bk]
 [7]  mg(t):mgrowth{exo}[log,bk]
-[8]  lam(t):productivity{exo}[log,bk]
-[9] epsg(t):prod_shock{iid}[]
-[10] epsu(t):mgro_shock{iid}[]
+[8]  epsg(t):prod_shock{iid}[]
+[9]  epsu(t):mgro_shock{iid}[]
 
 %Boundary Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 None
@@ -55,15 +54,6 @@ None
 [9]   @Fh(t+1)    = FF_1{@Fh(t)};
 [10]  @Fh_bar     = SS{@Fh(t)};
 
-# Utility function and derivatives
-[11]  @U(t)       = LOG(c(t))+B*h(t);
-[12]  @MU_c(t)    = DIFF{@U(t),c(t)};
-[13]  @MU_c_bar   = SS{@MU_c(t)};
-[14]  @MU_c(t+1)  = FF_1{@MU_c(t)};
-[15]  @MU_h(t)    = DIFF{@U(t),h(t)};
-[16]  @MU_h_bar   = SS{@MU_h(t)};
-[17]  @MU_h(t+1)  = FF_1{@MU_h(t)};
-
 
 
 %Non-Linear First-Order Conditions++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -72,17 +62,14 @@ None
 # The economy's budget constraint
 [1]  k(t)+(1/p(t))-(1-delta)*k(t-1)-w(t)*h(t)-r(t)*k(t-1) = 0;
 
-# Cash-in-advance constraint
-[2]  p(t)*c(t)-1 = 0;
-
 # Defining the real return on physical capital
 [3]  r(t)-@Fk(t) = 0;
 
 # Defining the economy's real wage
 [4]  w(t)-@Fh(t) = 0;
 
 # The MRS between consumption and leisure
-[5]  (B/(w(t)*p(t)))+(betta/(E(t)|p(t+1)*E(t)|c(t+1)*mg(t))) = 0;
+[5]  (B/(w(t)*p(t)))+(betta/mg(t)) = 0;
 
 # The asset equation for bonds
 [6]  betta*(w(t)/E(t)|w(t+1))*(1+E(t)|r(t+1)-delta)-1 = 0;
@@ -95,14 +82,12 @@ None
 
 
 %Steady States [Closed Form]+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-[2]  r_bar   = (1/betta)-(1-delta);
-[3]  w_bar   = (1-theta)*(r_bar/theta)**(theta/(theta-1.0));
-[4]  c_bar   = -(betta*w_bar)/(mg_bar*B);
-[5]  p_bar   = 1/c_bar;
-[6]  k_bar   = c_bar/((r_bar/theta)-delta);
-[7]  h_bar   = (r_bar/theta)**(1/(1-theta))*k_bar;
-[8]  l_bar   = 1.0-h_bar;
-[9]  y_bar   = c_bar+delta*k_bar;
+[1]  r_bar   = (1/betta)-(1-delta);
+[2]  w_bar   = (1-theta)*(r_bar/theta)**(theta/(theta-1.0));
+[3]  c_bar   = -(betta*w_bar)/(mg_bar*B);
+[4]  p_bar   = 1/c_bar;
+[5]  k_bar   = c_bar/((r_bar/theta)-delta);
+[6]  h_bar   = (r_bar/theta)**(1/(1-theta))*k_bar;
 
 
 %Steady State Non-Linear System [Manual]+++++++++++++++++++++++++++++++++++++++++++++++++