/pymaclab/dsge/macrolab.py

Large files files are truncated, but you can click here to view the full file

'''
.. module:: macrolab
   :platform: Linux
   :synopsis: The macrolab module contains the the most important classes for doing work with DSGE models. In particular it supplies
              the DSGEmodel class containing most of the functionality of DSGE model instances. The module also contains the TSDataBase
              class which was supposed to be an advanced data carrier class to be passed to the DSGE model instance for estimation
              purposes, but this is deprecated and will probably be replaced with a pandas data frame in the near future.

.. moduleauthor:: Eric M. Scheffel <eric.scheffel@nottingham.edu.cn>


'''

#from __future__ import division
import copy
from copy import deepcopy
import os
import sys
import re
from helpers import now_is
import numpy as np
from numpy import matlib as mat
from scipy.linalg import eig as scipyeig
from scipy import optimize
from tempfile import mkstemp
import time
import datetime
import glob
# Import Sympycore, now comes supplied and installed with pymaclab
import sympycore
# Import PP, now comes supplied and installed with pymaclab
import pp

#NOTE: Imports from the refactor
from pymaclab.filters._hpfilter import hpfilter
from pymaclab.filters._bkfilter import bkfilter
from ..stats.var import VAR
from solvers.steadystate import SSsolvers, ManualSteadyState, Fsolve
from parsers._modparser import parse_mod
from parsers._dsgeparser import populate_model_stage_one,populate_model_stage_one_a,\
     populate_model_stage_one_b,populate_model_stage_one_bb,populate_model_stage_two
from updaters.one_off import Updaters, dicwrap, dicwrap_deep, listwrap, matwrap
from updaters.queued import Updaters_Queued, dicwrap_queued, dicwrap_deep_queued, listwrap_queued,\
     matwrap_queued, Process_Queue, queue

#Define a list of the Greek Alphabet for Latex
greek_alph = ['alpha','beta','gamma','delta','epsilon',
          'varepsilon','zeta','eta','theta','vartheta',
          'iota','kappa','lambda','mu','nu','xi','pi',
          'varpi','rho','varrho','sigma','varsigma',
          'tau','upsilon','phi','varphi','chi','psi',
          'omega','Gamma','Delta','Theta','Lambda',
          'Xi','Pi','Sigma','Upsilon','Phi','Psi','Omega']

# Setting paths
root = ''
modfpath = ''
datapath = ''
browserpath = ''
txtedpath = ''
texedpath = ''
pyedpath = ''
pdfpath = ''
lapackpath = ''
lapackname = ''
mlabpath = ''
# Is matlab installed and can it therefore be used?
#TODO: this is not safe at all if matlab is not installed
try:
    import mlabraw
    use_matlab = True
except:
    use_matlab = False
    sess1 = None
# Open mlab session
if use_matlab:
    sess1 = mlabraw.open('matlab - nojvm -nosplash')
    
# Shell class for derivatives branch
class Derivatives(object):
    pass

# Shell class for derivatives branch
class Inits(object):
    
    def __init__(self,other=None):
        self._other = other

    def init1(self,no_wrap=False):
        '''
        The init1 method. Model population proceeds from the __init__function here. In particular the data gets read in
        (not implemented at the moment) and the model parsing begins.
        
        .. note::
           The other.vardic gets created and the manual as well as
           the numerical steady state sections gets parsed and attached to the DSGE model instance. So the most import fields created
           here using function populate_model_stage_one() are:
           
             * other.vardic - variable names dictionary with transform and filtering info
             * other.mod_name - short name of the DSGE model from mod file
             * other.mod_desc - longer model description from mod file
             * other.paramdic - dic of defined parameters with their numerical values
             * other.manss_sys - list of equations from the closed form steady state section
             * other.ssys_list - list of equations from the numerical steady state section
           
           Also the updaters and updaters_queued branches are opened here and the other.vardic gets wrapped for dynamic updating
           behaviour.
        
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst
        
        '''
        other = self._other
        initlev = other._initlev
        ncpus = other._ncpus
        mk_hessian = other._mk_hessian
        mesg = other._mesg
        # Attach the data from database
        if other.dbase != None:
            other.getdata(dbase=other.dbase)        
        if mesg: print "INIT: Instantiating DSGE model with INITLEV="+str(initlev)+" and NCPUS="+str(ncpus)+"..."
        if mesg: print "INIT: Attaching model properties to DSGE model instance..."

        # Create None tester regular expression
        #        _nreg = '^\s*None\s*$'
        #        nreg = re.compile(_nreg)
        
        if mesg: print "INIT: Parsing model file into DSGE model instance..."
        txtpars = parse_mod(other.modfile)
        other.txtpars = txtpars  # do we need txtpars attached for anything else?
        secs = txtpars.secs # do we need txtpars attached for anything else?
        if mesg: print "INIT: Extraction of info into DSGE model instance Stage [1]..."
        # Initial population method of model, does NOT need steady states
        other = populate_model_stage_one(other, secs)

        if not no_wrap:
            # Open updaters path
            other.updaters = Updaters()        
            # Open the updaters_queued path
            other.updaters_queued = Updaters_Queued()
            # Add the empty queue
            other.updaters_queued.queue = queue
    
            # Wrap the vardic
            other.updaters.vardic = dicwrap_deep(other,'self.vardic',initlev)
            other.updaters_queued.vardic = dicwrap_deep_queued(other,'self.vardic',initlev)
        
    def init1a(self):
        '''
        The init1a method. Model population proceeds from the init1 method here. The only field which get created here is the raw
        (i.e. unsubstituted) substitution dictionary.
        
        .. note::
            Field which are created here using the function populate_model_stage_one_a() which in turn calls mk_subs_dic():
           
             * other.nlsubs_raw1 - a list of the @items and their replacements
             * other.nlsubsdic - the above just expressed as a keyed list
        
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst

        '''
        other = self._other
        mesg = other._mesg
        secs = other.txtpars.secs
        # This is an additional populator which creates subsitution dictionary
        # and uses this to already get rid of @s in the manuall sstate list
        if mesg: print "INIT: Extraction of info into DSGE model instance Stage [2]..."
        # This function only creates the raw substitution dictionary and list from modfile
        other = populate_model_stage_one_a(other,secs)
        
    def init1b(self,no_wrap=False):
        '''
        The init1b method. Model population proceeds from the init1a method here.

        .. note::
        
           The only thing which gets done here purposefully *after* calling init1a() is to wrap these fields:
           
             * other.nlsubsdic - the above just expressed as a keyed list
             * other.paramdic - dic of defined parameters with their numerical values
             
           in order to give them dynamic updating behaviour. No more is done in this init method call.
        
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst

        '''
        other = self._other
        initlev = other._initlev
        secs = other.txtpars.secs
        other = populate_model_stage_one_b(other,secs)
        
        if not no_wrap:
            # Wrap the nlsubsdic
            if 'nlsubsdic' in dir(other): other.updaters_queued.nlsubsdic = dicwrap_queued(other,'self.nlsubsdic',initlev)
            # Wrap the paramdic
            other.updaters_queued.paramdic = dicwrap_queued(other,'self.paramdic',initlev)
            
    
            # Wrap the nlsubsdic
            if 'nlsubsdic' in dir(other): other.updaters.nlsubsdic = dicwrap(other,'self.nlsubsdic',initlev)
            # Wrap the paramdic
            other.updaters.paramdic = dicwrap(other,'self.paramdic',initlev)        
        
    def init1c(self,no_wrap=False):
        '''
        The init1c method. Model population proceeds from the init1b method here. We call populate_model_stage_one_bb() which does quite
        a bit of substitution/replacement of the @-prefixed variables. It does this in the numerical and closed form steady state
        calculation sections, as well as in the actual FOCs themselves.

        .. note::
        
           *After* that we can wrap various fields for updating which are:
           
             * other.foceqs - list of firs-order conditions with @ replacements done
             * other.manss_sys - the list of equations from closed form SS section
             * other.syss_list - the list of equations from numerical SS section
        
           Notice also that here we replace or generate fields in case the FOCs are supposed to be used
           directly in SS calculation because the "use_focs" parameter was not passed empty.
        
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst

        '''
        other = self._other
        initlev = other._initlev
        secs = other.txtpars.secs
        mesg = other._mesg
        if mesg: print "INIT: Substituting out @ variables in steady state sections..."
        other = populate_model_stage_one_bb(other,secs)
        
        if not no_wrap:
            if 'foceqs' in dir(other):
                # Wrap foceqs
                other.updaters.foceqs = listwrap(other,'self.foceqs',initlev)
                # Wrap foceqs
                other.updaters_queued.foceqs = listwrap(other,'self.foceqs',initlev)
        
        # Allow for numerical SS to be calculated using only FOCs
        if other._use_focs and other._ssidic != None:
            list_tmp = []
            for elem in other._use_focs:
                list_tmp.append(other.foceqss[elem])
            other.ssys_list = deepcopy(list_tmp)
            # Check if this has not already been produced in populate_model_stage_one_bb, chances are it has
            if 'ssidic' not in dir(other): other.ssidic = copy.deepcopy(other._ssidic)
        
        if not no_wrap:   
            # Wrap manss_sys
            if 'manss_sys' in dir(other):
                other.updaters.manss_sys = listwrap(other,'self.manss_sys',initlev)
                other.updaters_queued.manss_sys = listwrap_queued(other,'self.manss_sys',initlev)
            # Wrap ssys_list
            if 'ssys_list' in dir(other):
                other.updaters.ssys_list = listwrap(other,'self.ssys_list',initlev)
                other.updaters_queued.ssys_list = listwrap_queued(other,'self.ssys_list',initlev)


    def init2(self):
        '''
        The init2 method. Model population proceeds from the init1c method here. In this initialisation method
        the only thing which is being done is to open up the sssolvers branch and pass down required objects to
        the manuals closed from solver and the numerical root-finding solver depending on whether information for
        this has been included in the DSGE model file. *No* attempt is made at solving for the steady state, the
        respective solvers are only being *prepared*.
        
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst
        
        '''
        other = self._other
        mesg = other._mesg
        secs = other.txtpars.secs
        initlev = other._initlev
        if mesg: print "INIT: Preparing DSGE model instance for steady state solution..."

        # Attach the steady state class branch, and add Fsolve if required but do no more !
        other.sssolvers = SSsolvers()

        # check if the Steady-State Non-Linear system .mod section has an entry
        if all([False if 'None' in x else True for x in secs['manualss'][0]]) or other._use_focs:
            intup = (other.ssys_list,other.ssidic,other.paramdic)
            other.sssolvers.fsolve = Fsolve(intup,other=other)
        # check if the Steady-State Non-Linear system closed form has an entry
        #  we do this here with ELIF because we only want to set this up for solving if manualss does not exist
        elif all([False if 'None' in x else True for x in secs['closedformss'][0]]):
            alldic = {}
            alldic.update(other.paramdic)
            intup = (other.manss_sys,alldic)
            other.sssolvers.manss = ManualSteadyState(intup,other=other)

    def init3(self):
        '''
        Initialisation method init3 is quite complex and goes through a number of logical tests in order to determine
        how to solve for the steady state of the model.
        
        .. note::
        
           There are 7 different ways a DSGE model can obtain its steady state solution depending on what information has
           been provided:
           
             1) The steady state values dictionary has been passed as argument, then init3() will NEVER be called
             2) Information has been provided using the "use_focs" parameter to use FOCs directly, externally passed using use_focs
             3) Information has been provided using the "use_focs" parameter to use FOCs directly, but inside model file
             4) Information has only been provided in the numerical SS section
             5) Information has only been provided in the closed form SS section
             6) Both CF-SS and NUM-SS info are present and NUM-SS is subset if CF-SS
             7) Both CF-SS and NUM-SS info are present and CF is residual
             
           These options are better explained in the documentation to PyMacLab in the steady state solver section.
           
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst
        
        '''
        other = self._other
        txtpars = other.txtpars
        secs = txtpars.secs
        initlev = other._initlev
################################## STEADY STATE CALCULATIONS !!! #######################################
        if other._mesg: print "INIT: Attempting to find DSGE model's steady state automatically..."
        # ONLY NOW try to solve !
        ##### OPTION 1: There is only information externally provided and we are using FOCs
        if other._use_focs and other._ssidic != None:
            if '_internal_focs_used' not in dir(other):
                if other._mesg: print "SS: Using FOCs and EXTERNALLY supplied information...attempting to solve SS..."
            else:
                ##### OPTION 1b: There is only information internally provided and we are using FOCs
                if other._mesg: print "SS: Using FOCs and INTERNALLY supplied information...attempting to solve SS..."
                del other._internal_focs_used
            other.sssolvers.fsolve.solve()
            if other.sssolvers.fsolve.ier == 1:
                other.sstate = deepcopy(other.sssolvers.fsolve.fsout)
                other.switches['ss_suc'] = ['1','1']
                if other._mesg: print "INIT: Steady State of DSGE model found (SUCCESS)..."
            else:
                other.switches['ss_suc'] = ['1','0']
                if other._mesg: print "INIT: Steady State of DSGE model not found (FAILURE)..."            
            return
            
        ##### OPTION 2: There is only information provided in the numerical section NOT in closed-form
        if not all([False if 'None' in x else True for x in secs['closedformss'][0]]) and\
           all([False if 'None' in x else True for x in secs['manualss'][0]]) and not other._use_focs and other._ssidic == None:
            if other._mesg: print "SS: ONLY numerical steady state information supplied...attempting to solve SS..."
            other.sssolvers.fsolve.solve()

        ##### OPTION 3: There is only information on closed-form steady state, BUT NO info on numerical steady state
        # Solve using purely closed form solution if no other info on model is available
        if all([False if 'None' in x else True for x in secs['closedformss'][0]]) and\
           not all([False if 'None' in x else True for x in secs['manualss'][0]]):
            if other._mesg: print "SS: ONLY CF-SS information supplied...attempting to solve SS..."
            alldic = {}
            alldic.update(other.paramdic)
            intup = (other.manss_sys,alldic)
            other.sssolvers.manss = ManualSteadyState(intup)
            other.sssolvers.manss.solve()
            other.sstate = {}
            other.sstate.update(other.sssolvers.manss.sstate)
        ##### OPTION 4: There is information on closed-form AND on numerical steady state
        # Check if the numerical and closed form sections have entries
        if all([False if 'None' in x else True for x in secs['manualss'][0]]) and\
           all([False if 'None' in x else True for x in secs['closedformss'][0]]):
            # Create unordered Set of closed from solution variables
            manss_set = set()
            for elem in other.manss_sys:
                manss_set.add(elem.split('=')[0].lstrip().rstrip())
            # If ssidic is not empty we need to make sure it perfectly overlaps with manss_set in order to replace ssidic
            if other.ssidic != {}:
                numss_set = set()
                for keyo in other.ssidic.keys():
                    numss_set.add(keyo)
                ##### OPTION 4a: If there is an ssidic and its keys are subset of manss_set, the use as suggestion for new ssi_dic
                if numss_set.issubset(manss_set):
                    if other._mesg: print "SS: CF-SS and NUM-SS (overlapping) information information supplied...attempting to solve SS..."
                    alldic = {}
                    alldic.update(other.paramdic)
                    intup = (other.manss_sys,alldic)
                    other.sssolvers.manss = ManualSteadyState(intup)
                    other.sssolvers.manss.solve()
                    for keyo in other.ssidic.keys():
                        other.ssidic[keyo] = other.sssolvers.manss.sstate[keyo]
            ##### OPTION 4b: ssidic is empty, so we have to assumed that the variables in closed form are suggestions for ssidic
            # If it is empty, then just compute the closed form SS and pass to ssidic as starting value
            elif other.ssidic == {}:
                if other._mesg: print "SS: CF-SS and NUM-SS (empty ssdic) information information supplied...attempting to solve SS..."
                alldic = {}
                alldic.update(other.paramdic)
                intup = (other.manss_sys,alldic)
                other.sssolvers.manss = ManualSteadyState(intup)
                other.sssolvers.manss.solve()
                other.ssidic.update(other.sssolvers.manss.sstate)
                # Test at least if the number of ssidic vars equals number of equations
                if len(other.ssidic.keys()) != len(other.ssys_list):
                    print "Error: Number of variables in initial starting values dictionary != number of equations"
                    sys.exit()
        ######## Finally start the numerical root finder with old or new ssidic from above
        if all([False if 'None' in x else True for x in secs['manualss'][0]]) and\
           not all([False if 'None' in x else True for x in secs['closedformss'][0]]) and not other._use_focs:            
            other.sssolvers.fsolve.solve()
        elif all([False if 'None' in x else True for x in secs['manualss'][0]]) and\
           all([False if 'None' in x else True for x in secs['closedformss'][0]]) and not other._use_focs:           
            other.sssolvers.fsolve.solve()

        if all([False if 'None' in x else True for x in secs['manualss'][0]]):
            if other.sssolvers.fsolve.ier == 1:
                other.sstate = other.sssolvers.fsolve.fsout
                other.numssdic = other.sssolvers.fsolve.fsout
                # Attach solutions to intial variable dictionaries, for further analysis
                other.ssidic_modfile = deepcopy(other.ssidic)
                # Update old ssidic with found solutions
                other.ssidic = deepcopy(other.sssolvers.fsolve.fsout)
                other.sssolvers.fsolve.ssi = other.ssidic
                other.switches['ss_suc'] = ['1','1']
                if other._mesg: print "INIT: Steady State of DSGE model found (SUCCESS)..."
            else:
                other.switches['ss_suc'] = ['1','0']
                if other._mesg: print "INIT: Steady State of DSGE model not found (FAILURE)..."

        ########## Here we are trying to merge numerical SS solver's results with result closed-form calculations, if required
        if all([False if 'None' in x else True for x in secs['manualss'][0]]) and\
           all([False if 'None' in x else True for x in secs['closedformss'][0]]):
            if other._mesg: print "INIT: Merging numerical with closed form steady state if needed..."
        # Check for existence of closedform AND numerical steady state
        # We need to stop the model instantiation IFF numerical solver was attempted but failed AND closed form solver depends on it.
        if all([False if 'None' in x else True for x in secs['closedformss'][0]]) and\
           all([False if 'None' in x else True for x in secs['manualss'][0]]) and other.ssidic != {}:
            if other.switches['ss_suc'] == ['1','0']:
                print "ERROR: You probably want to use numerical steady state solution to solve for RESIDUAL closed form steady states."
                print "However, the numerical steady state solver FAILED to find a root, so I am stopping model instantiation here."
                sys.exit()
            ##### OPTION 5: We have both numerical and (residual) closed form information
            # Check if a numerical SS solution has been attempted and succeeded, then take solutions in here for closed form.
            elif other.switches['ss_suc'] == ['1','1'] and not numss_set.issubset(manss_set):
                if other._mesg: print "SS: CF-SS (residual) and NUM-SS information information supplied...attempting to solve SS..."
                alldic = {}
                alldic.update(other.sstate)
                alldic.update(other.paramdic)
                intup = (other.manss_sys,alldic)
                other.sssolvers.manss = ManualSteadyState(intup)
                other.sssolvers.manss.solve()
                # Here merging takes place
                other.sstate.update(other.sssolvers.manss.sstate)

        # Double check if no steady state values are negative, as LOGS may have to be taken.
        if 'sstate' in dir(other):
            for keyo in other.sstate.keys():
                if '_bar' in keyo and float(other.sstate[keyo]) < 0.0:
                    print "WARNING: Steady state value "+keyo+ " is NEGATIVE!"
                    print "This is very likely going to either error out or produce strange results"
                    print "Re-check your model declarations carefully!"
######################################### STEADY STATE CALCULATION SECTION DONE ##################################
    def init4(self,no_wrap=False):
        '''
        This model instance sub-initializor only calls the section which use the computed steady state
        in order to compute derivatives and open dynamic solver branches on the instance. But Jacobian and Hessian
        are *not* computed here, this is postponed to the next init level.
        
        .. note::
        
           The following last field is wrapped for dynamic execution:
           
             * other.sigma - the variance-covariance matrix of the iid shocks
             
           Notice that after wrapping this last field the process_queue class is instantiated at last,
           because it needs to have access to *all* of the wrapped fields. Also in this method, the function
           populate_model_stage_two() is called which prepares the nonlinear FOCs for derivative-taking.
           
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst

        '''
        other = self._other
        txtpars = other.txtpars
        secs = txtpars.secs
        initlev = other._initlev
        ncpus = other._ncpus
        mk_hessian = other._mk_hessian
        mesg = other._mesg
        if mesg: print "INIT: Preparing DSGE model instance for computation of Jacobian and Hessian..."
        # Now populate more with stuff that needs steady state
        other = populate_model_stage_two(other, secs)

        if not no_wrap:
            # Need to wrap variance covariance matrix here
            other.updaters.sigma = matwrap(other,'self.sigma',initlev)
            # Need to wrap variance covariance matrix here
            other.updaters_queued.sigma = matwrap_queued(other,'self.sigma',initlev)

            ####### All queued updaters initialized, no add processing instance
            # Add the queue process instance
            other.updaters_queued.process_queue = Process_Queue(other=other)

    def init5(self,update=False):
        '''
        This model instance initialisation step is the last substantial one in which the dynamic solution of the DSGE
        model instance is finally computed using a choice of methods which can be called at runtime.
        
         
        :param other:     The DSGE model instance itother.
        :type other:      dsge_inst

        '''
        other = self._other
        other.derivatives = Derivatives()
        txtpars = other.txtpars
        secs = txtpars.secs
        initlev = other._initlev
        ncpus = other._ncpus
        mk_hessian = other._mk_hessian
        mesg = other._mesg
        #TODO: delay above and only import if needed
        from solvers.modsolvers import MODsolvers
        # Open the model solution tree branch
        other.modsolvers = MODsolvers()
        ######################## LINEAR METHODS !!! ############################
        # see if there are any log-linearized equations
        if all([False if 'None' in x else True for x in secs['modeq'][0]]):
            from solvers.modsolvers import PyUhlig, MatUhlig, MatKlein, MatKleinD, ForKlein
            if mesg: print "INIT: Computing DSGE model's log-linearized solution using Uhlig's Toolbox..."

            # Open the matlab Uhlig object
            intup = ((other.nendo,other.ncon,other.nexo),
                 other.eqindx,
                 other.vreg,
                 other.llsys_list,
                 other.diffli1,
                 other.diffli2,
                 sess1,
                 other.vardic)
            other.modsolvers.matuhlig = MatUhlig(intup)

            # Open the native Uhlig object
            intup = ((other.nendo,other.ncon,other.nexo),
                 other.eqindx,
                 other.vreg,
                 other.llsys_list,
                 other.diffli1,
                 other.diffli2,
                 sess1)
            other.modsolvers.pyuhlig = PyUhlig(intup)

            # Open the matlab Klein object
            intup = ((other.nendo,other.ncon,other.nexo),
                 other.eqindx,
                 other.vreg,
                 other.llsys_list,
                 other.diffli1,
                 other.diffli2,
                 sess1)
            other.modsolvers.matklein = MatKlein(intup)

            # Open the Fortran Klein object
            intup = ((other.nendo,other.ncon,other.nexo),
                 other.eqindx,
                 other.vreg,
                 other.llsys_list,
                 other.diffli1,
                 other.diffli2,
                 sess1)
            other.modsolvers.forklein = ForKlein(intup)
    ################## 1ST-ORDER NON-LINEAR METHODS !!! ##################
        if all([False if 'None' in x else True for x in secs['focs'][0]]):
            from solvers.modsolvers import (MatWood, ForKleinD)
            
            if not update:
                if ncpus > 1 and mk_hessian:
                    if mesg: print "INIT: Computing DSGE model's Jacobian and Hessian using parallel approach..."
                    other.mkjahepp()
                elif ncpus > 1 and not mk_hessian:
                    if mesg: print "INIT: Computing DSGE model's Jacobian using parallel approach..."
                    other.mkjahepp()
                else:
                    if mesg: print "INIT: Computing DSGE model's Jacobian and Hessian using serial approach..."
                    other.mkjahe()
            else:
                if ncpus > 1 and mk_hessian:
                    if mesg: print "INIT: Computing DSGE model's Jacobian and Hessian using parallel approach..."
                    other.mkjaheppn()
                elif ncpus > 1 and not mk_hessian:
                    if mesg: print "INIT: Computing DSGE model's Jacobian using parallel approach..."
                    other.mkjaheppn()
                else:
                    if mesg: print "INIT: Computing DSGE model's Jacobian and Hessian using serial approach..."
                    other.mkjahen()                

            # Check if the obtained matrices A and B have correct dimensions
            if other.derivatives.jAA.shape[0] != other.derivatives.jAA.shape[1]:
                print "ERROR: Matrix A of derivatives does not have #vars=#equations"
            if other.derivatives.jBB.shape[0] != other.derivatives.jBB.shape[1]:
                print "ERROR: Matrix B of derivatives does not have #vars=#equations"
            # Open the MatWood object
            intup = (other.derivatives.jAA,other.derivatives.jBB,
                 other.nexo,other.ncon,
                 other.nendo,sess1)
            other.modsolvers.matwood = MatWood(intup)
            # Make the AA and BB matrices as references available instead
            other.modsolvers.matwood.jAA = other.derivatives.jAA
            other.modsolvers.matwood.jBB = other.derivatives.jBB
            # Open the Fortran KleinD object
            if 'nlsubsys' in dir(other):
                intup = (other.derivatives.numj,
                     other.nendo,other.nexo,
                     other.ncon,other.sigma,
                     other.derivatives.jAA,other.derivatives.jBB,
                     other.vardic,other.vdic,
                     other.mod_name,other.audic,
                     other.derivatives.numjl,
                     other.nother)
            else:
                intup = (other.derivatives.numj,
                     other.nendo,other.nexo,
                     other.ncon,other.sigma,
                     other.derivatives.jAA,other.derivatives.jBB,
                     other.vardic,other.vdic,
                     other.mod_name,other.audic)
            other.modsolvers.forkleind = ForKleinD(intup,other=other)
            # Make the AA and BB matrices as references available instead
            other.modsolvers.forkleind.A = other.derivatives.jAA
            other.modsolvers.forkleind.B = other.derivatives.jBB

    ################## 2ND-ORDER NON-LINEAR METHODS !!! ##################
        if all([False if 'None' in x else True for x in secs['vcvm'][0]]) and 'numh' in dir(other.derivatives):
            from solvers.modsolvers import (PyKlein2D, MatKlein2D)
            # Open the MatKlein2D object
            if 'nlsubsys' in dir(other):
                intup = (other.derivatives.numj,other.derivatives.numh,
                     other.nendo,other.nexo,
                     other.ncon,other.sigma,
                     other.derivatives.jAA,other.derivatives.jBB,
                     other.vardic,other.vdic,
                     other.mod_name,other.audic,
                     other.derivatives.numjl,other.derivatives.numhl,
                     other.nother,sess1)
            else:
                intup = (other.derivatives.numj,other.derivatives.numh,
                     other.nendo,other.nexo,
                     other.ncon,other.sigma,
                     other.derivatives.jAA,other.derivatives.jBB,
                     other.vardic,other.vdic,
                     other.mod_name,other.audic,
                     sess1)
            other.modsolvers.matklein2d = MatKlein2D(intup)
            # Make the AA and BB matrices as references available instead
            other.modsolvers.matklein2d.A = other.derivatives.jAA
            other.modsolvers.matklein2d.B = other.derivatives.jBB

            # Open the PyKlein2D object, but don't pass mlabwrap session
            intup = intup[:-1]
            other.modsolvers.pyklein2d = PyKlein2D(intup,other=other)
            # Make the AA and BB matrices as references available instead
            other.modsolvers.pyklein2d.A = other.derivatives.jAA
            other.modsolvers.pyklein2d.B = other.derivatives.jBB
            other.modsolvers.pyklein2d.forkleind.A = other.derivatives.jAA
            other.modsolvers.pyklein2d.forkleind.B = other.derivatives.jBB
            
    def init_out(self):
        '''
        The final intializor section does some extra stuff after all has been done.
           
        :param self:     The DSGE model instance itself.
        :type self:      dsge_inst

        '''
        other = self._other
        initlev = other._initlev
        # Make sure the jobserver has done his jobs
        jobserver.wait()

# Empty locdic for the locate helper function
locdic = {}

##################THE DSGE MODEL CLASS (WORKS)#####################
class DSGEmodel(object):
    '''
    This is the macrolab DSGEmodel class. It is the main class of the packages and instantiates
    DSGE model instances which possess many features such as model file parsers, solvers, etc. The __init__
    function first called mostly attaches the passed arguments to the DSGE model instance in form of private _X
    data fields.
    
    .. note::
    
       Notice that the various init method, i.e. init1, init1a, etc. are not called here, but they are called externally in the
       pymaclab package when instantiating a new DSGE model using pymaclab.newMOD().
    
    :param ffile:             The absolute path to the PyMacLab DSGE model file to be parsed on instantiation
    :type ffile:              str
    :param dbase:             A database with time series data for estimation purposes. Not implemented at the moment
    :type standard:           unknown
    :param initlev:           Takes values 0,1,2. Determines how deep the instantiation cascades through all methods.
                              If 0 then the file is parsed and the instance is only *prepared* for steady state solving.
                              If 1 then the file is parsed, the SS is automatically computed and instance is *prepared*
                              for dynamic solving. If 2 then the instance is solved all the way.
    :type initlev:            int
    :param mesg:              If True then lots of diagnostics are printed to the screen during instantiation.
    :type mesg:               bool
    :param ncpus:             The number of CPU core to be employed, defaults to 1. But 'auto' can also be used for detection
    :type ncpus:              int|str
    :param mk_hessian:        Whether the Hessian should be computed as this is expensive.
    :type mk_hessian:         bool
    :param use_focs:          Should the FOCs be used directly to look for the steady state? Must use list|tuple to pick equations.
    :type use_focs:           tuple
    :param ssidic:            A Python ssidic with the initial starting values for solving for the SS numerically
    :type ssidic:             dic
    :param sstate:            A Python ssidic with the externally computed steady state values of the model
    :type sstate:             dic
        
    :return self:             *(dsge_inst)* - A populated DSGE model instance with fields and methods

    '''
    # Instantiate a global class jobserver accessible from all instances
    global ppservers, jobserver
    ppservers = ()
    jobserver = pp.Server(ppservers=ppservers)
    # Initializes the absolute basics, errors difficult to occur
    def __init__(self,ffile=None,dbase=None,initlev=2,mesg=False,ncpus='auto',mk_hessian=True,\
                 use_focs=False,ssidic=None,sstate=None,vtiming={'exo':[-1,0],'endo':[-1,0],'con':[0,1]}):
        # Open the inits branch
        self.inits = Inits(other=self)
        # Make sure the jobserver has done his global jobs each time you instantiate a new instance
        jobserver.wait()
        if sstate != None:
            self._sstate = deepcopy(sstate)
        self._initlev = initlev #TODO: remove this as an option
        self._vtiming = vtiming
        self._mesg = mesg
        self._ncpus = ncpus
        self._mk_hessian = mk_hessian
        self._use_focs = use_focs
        self._ssidic = ssidic
        if self._use_focs and self._ssidic != None:
            self.ssidic = deepcopy(self._ssidic)
        # Set no author
        self.setauthor()
        # Open the switches dic and initiate
        #NOTE this is _not_ pythonic
        self.switches = {}
        self.switches['errocc'] = '0'
        self.switches['ss_suc'] = ['0','0']
        # Attach some model attributes
        self.modfile = ffile
        self.dbase = dbase

        
    def mk_dynare(self,order=1,centralize=False,fpath=None,focli=None):
        # Import the template and other stuff
        from pymaclab.modfiles.templates import mako_dynare_template
        from scipy import io
        from solvers.modsolvers import Dynare
        import tempfile
        import os
        import glob
        template_dic = deepcopy(self.template_paramdic)
        
        # Render the template to be passed to dynare++
        tmp_dic = {}
        tmp_dic['focli'] = focli
        template_dic.update(tmp_dic)
        modstr = mako_dynare_template.render(**template_dic)
        
        # If a filepath has been passed then just write the Dynare++ modfile, but no more!
        if fpath != None:
            filo = open(fpath,'w')
            filo.write(modstr)
            filo.flush()
            filo.close()
            return

        filo = tempfile.NamedTemporaryFile()
        fname = filo.name
        fname2 = fname.split('/')[-1]
        filo2 = open(os.path.join(os.getcwd(),'test.mod'),'w')
        filo2.write(modstr)
        filo2.flush()
        filo2.close()
        filo.file.write(modstr)
        filo.file.flush()
        self.modsolvers.dynare = Dynare({})
        self.modsolvers.dynare.__setattr__('modfile',modstr)
        if not centralize:
            os.system('dynare++ --no-centralize --order '+str(order)+' '+filo.name)
        else:
            os.system('dynare++ --order '+str(order)+' '+filo.name)
        dynret = io.loadmat(os.path.join(os.getcwd(),filo.name.split('/')[-1]+'.mat'))
        # Check if solution has been computed and attache all solution matrices to dynare instance
        if dynret.has_key('dyn_g_1'):
            self.modsolvers.dynare.__init__(dynret)
        else:
            print "FAIL: Dynare could not determine solution."
        filo.close()
        # Delete all of the files
        for filor in glob.glob(fname2+'*.*'):
            os.remove(filor)

        # Create some other objects on dynare branch
        self.modsolvers.dynare.sstate = {}
        self.modsolvers.dynare.dynorder = []
        self.modsolvers.dynare.dynsorder = []
        for i1,elem in enumerate(self.modsolvers.dynare.dyn_vars):
            self.modsolvers.dynare.dynorder.append(str(elem).strip()+'(t)')
            if str(elem).strip()+'_bar' in self.sstate.keys():
                self.modsolvers.dynare.sstate[str(elem).strip()+'_bar'] = self.modsolvers.dynare.dyn_ss[i1,0]
        for i1,elem in enumerate(self.modsolvers.dynare.dyn_state_vars):
            self.modsolvers.dynare.dynsorder.append(str(elem).strip()+'(t)')

        # Create P and F matrix reflecting ordering of pymaclab
        P = np.matlib.zeros([self.nstat,self.nstat])
        F = np.matlib.zeros([self.ncon,self.nstat])
        listo = [x[0] for x in self.vdic['exo']]+[x[0] for x in self.vdic['endo']]
        for i1,elem in enumerate(listo):
            for i2,elem2 in enumerate(listo):
                P[i1,i2] = self.modsolvers.dynare.dyn_g_1[self.modsolvers.dynare.dynorder.index(elem),self.modsolvers.dynare.dynsorder.index(elem2)]
        for i1,elem in enumerate([x[0] for x in self.vdic['con']]):
            for i2,elem2 in enumerate(listo):
                F[i1,i2] = self.modsolvers.dynare.dyn_g_1[self.modsolvers.dynare.dynorder.index(elem),self.modsolvers.dynare.dynsorder.index(elem2)]
        self.modsolvers.dynare.PP = deepcopy(P)
        self.modsolvers.dynare.FF = deepcopy(F)
            
    def find_rss(self,mesg=False,rootm='hybr',scale=0.0):
        '''
        The is a method which can be called to find the risky steady state
           
        :param self:     The DSGE model instance itself.
        :type self:      dsge_inst

        '''
         
        varbar = []
        nexo = self.nexo
        for elem in self.vardic['endo']['var']:
            varbar.append(elem[0].split('(')[0]+'_bar')
        for elem in self.vardic['con']['var']:
            varbar.append(elem[0].split('(')[0]+'_bar')
        tmp_dic = {}
        tmp_dic.update(self.paramdic)
        tmp_dic.update(self.sstate)
        sstate_li = []
        sstate = {}
        for elem in varbar:
            if elem in tmp_dic.keys():
                sstate[elem] = tmp_dic[elem]
                sstate_li.append(tmp_dic[elem])
        rsstate_li = deepcopy(sstate_li)
        rsstate = deepcopy(sstate)

        clone = copy.deepcopy(self)
        clone._mesg = False
        
        # Define the function to be handed over
        # to fsolve
        def func(invar):
            '''
            ####### Put in a trap in case number turn negative ########
            negli = [True if x < 0.0 else False for x in invar]
            if any(negli):
                for i1,elem in enumerate(negli):
                    invar[i1] = self.sstate[varbar[i1]]+0.01*self.sstate[varbar[i1]]
            ###########################################################
            '''
            # Update ss with passed values
            for i1,elem in enumerate(varbar):            
                clone.sstate[elem] = invar[i1]
            # Update the model's derivatives, but only numerically
            clone.init5(update=True)
            # Solve the 2nd-order accurate solution
            clone.modsolvers.pyklein2d.solve()
            # Get retval into right shape
            KX = clone.modsolvers.pyklein2d.KX
            retval = [float(x) for x in KX[nexo:]]
            KY = clone.modsolvers.pyklein2d.KY
            retval+=[float(x) for x in KY]
            retval = np.array(retval)
            if mesg:
                print invar
                print '------------------------------------'
            return retval

        
        # Define the initial values and
        # start the non-linear solver
        init_val = np.array([float(self.sstate[varbar[i1]])+\
                             (scale/100.0)*float(self.sstate[varbar[i1]]) for i1 in range(len(varbar))])
        outobj = optimize.root(func,init_val,method=rootm)
        rss_funcval1 = outobj.fun
        output = outobj.x
        mesg = outobj.message
        ier = outobj.status
        
        self.rss_funcval1 = rss_funcval1
                
        # Check the difference and do bounded minimisation (root-finding)
        diff_dic = {}
        bounds_dic = {}
        for i1,keyo in enumerate(varbar):
            if keyo in self.sstate.keys():
                if output[i1] != self.sstate[keyo]:
                    diff_dic[keyo] = output[i1]
                    bounds_dic[keyo] = (output[i1],output[i1])
        if bounds_dic != {}:
            # Also add the non-bar SS values to bounds_dic, before doing constrained minimisation
            for keyo in self.sstate.keys():
                if "_bar" not in keyo: bounds_dic[keyo] = (self.sstate[keyo],self.sstate[keyo])
            # Do constrained root finding here
            if 'fsolve' in dir(clone.sssolvers):
                clone.sssolvers.fsolve.solve(bounds_dic=bounds_dic)
            # Also need to make sure residually computed SS get updated accordingly
            if 'manss' in dir(clone.sssolvers) and 'fsolve' in dir(clone.sssolvers):
                    clone.sssolvers.manss.paramdic.update(clone.sssolvers.fsolve.fsout)
                    clone.sssolvers.manss.solve()

            # Run loop one last time to get rss_funval2
            if 'fsolve' in dir(clone.sssolvers):
                clone.sstate.update(clone.sssolvers.fsolve.fsout)
            if 'manss' in dir(clone.sssolvers):
                clone.sstate.update(clone.sssolvers.manss.sstate)
            clone.init5(update=True)
            clone.modsolvers.pyklein2d.solve()
            retval = []
            # Get retval into right shape
            KX = clone.modsolvers.pyklein2d.KX
            retval = [float(x) for x in KX]
            KY = clone.modsolvers.pyklein2d.KY
            for elem in KY:
                retval.append(float(elem))
            retval = np.array(retval)
            self.rss_funcval2 = retval
             
        # Now attach final results to instance
        self.rsstate = deepcopy(self.sstate)
        self.rparamdic = deepcopy(self.paramdic)
        for i1,elem in enumerate(varbar):
            if elem in self.rsstate.keys(): self.rsstate[elem] = output[i1]
            if elem in self.rparamdic.keys(): self.rparamdic[elem] = output[i1]
        if bounds_dic != {}:
            if 'fsolve' in dir(clone.sssolvers):
                self.rsstate.update(clone.sssolvers.fsolve.fsout)
            if 'fsolve' in dir(clone.sssolvers) and 'manss' in dir(clone.sssolvers):
                self.rsstate.update(clone.sssolvers.manss.sstate)

            
        

    # html info of model opened with webbrowser
    def info(self):
        '''
        A convience method for collecting and displaying the model's properties in a browser
        using html language.
        '''
        tmplist = glob.glob('tempz*.html')
        for x in tmplist:
            os.remove(x)
        modname = self.mod_name
        secs = self.txtpars.secs
        direc = os.getcwd()
        fd,fpath = mkstemp(prefix='tempz',suffix='.html',dir=direc)
        htmlfile = os.fdopen(fd,'w+b')
        htmlfile.write('<HTML><BODY BGCOLOR="white">\n')
        htmlfile.write('<H2>%s</H2>\n'%modname)
        htmlfile.write('\n\n')
        for x1 in secs['info'][1]:
            htmlfile.write('<P>'+x1+'\n')
        htmlfile.write('<P>'+'<H4>Model Parameters</H4>\n')
        for x1 in secs['params'][1]:
            htmlfile.write('<P>'+x1+'\n')
        htmlfile.write('<P>'+'<H4>Nonlinear First-Order Conditions</H4>\n')
        for x1 in secs['focs'][1]:
            if x1[0] == '[':
                htmlfile.write('<P>'+x1+'\n')
            elif x1[0] != '[':
                htmlfile.write('<P>'+4*'&nbsp;'+ x1+'\n')
        htmlfile.write('<P>'+'<H4>Log-Linearized Model Equations</H4>\n')
        for x1 in secs['modeq'][1]:
            if x1[0] == '[':
                htmlfile.write('<P>'+x1+'\n')
            elif x1[0] != '[':
                htmlfile.write('<P>'+4*'&nbsp;'+ x1+'\n')
        htmlfile.write('</BODY></HTML>\n')
        htmlfile.write('<P>'+'<H2>Solution Section</H2>\n')
        htmlfile.write('<P>'+'<H4>Uhlig Toolkit Output</H4>\n')
        htmlfile.close()
        cmd = browserpath+' '+fpath
        os.system(cmd)
        tmplist = glob.glob('tempz*.html')
        for x in tmplist:
            os.remove(x)
        return 'Model Website opened!'
    
    def setauthor(self,author=None):
        '''
        Convience method for setting the model's author's name.
        '''
        if author == None:
            self.author = 'No author'
        else:
            self.author = author

    def pdf(self):
        '''
        This will pdflatex the model's tex file and then view it in the system's pdf viewer.
        '''
        modfile = self.modfile
        modfname = modfile.split('.')[0]
        if modfname+'.pdf' in os.listdir(modfpath):
            os.remove(os.path.join(modfpath,modfname+'.pdf'))
        if modfname+'.log' in os.listdir(modfpath):
            os.remove(os.path.join(modfpath,modfname+'.log'))
        if modfname+'.aux' in os.listdir(modfpath):
            os.remove(os.path.join(modfpath,modfname+'.aux'))
        if modfname+'.dvi' in os.listdir(modfpath):
            os.remove(os.path.join(modfpath,modfname+'.dvi'))
        if modfname+'.log' in os.listdir(os.getcwd()):
            os.remove(os.path.join(os.getcwd(),modfname+'.log'))

        # Does the model tex file even exist?
        if modfname+'.tex' not in os.listdir(modfpath):
            print 'Error: The model tex file does not exist!'
            print 'Use model.texed() to create a new one!'
            return

        # First check for Chktex syntax errors
        if 'texer.log' in os.listdir(os.getcwd()):
            os.remove(os.path.join(os.getcwd(),'texer.log'))
        args = '--quiet -v2 -n3 -n25 -n12 -n35'
        cmd = 'chktex '+args+' '+os.path.join(modfpath,modfname+'.tex'+' > texer.log')
        os.system(cmd)
        file = open(os.path.join(os.getcwd(),'texer.log'),'rU')
        errlist = file.readlines()
        if len(errlist)== 0:
            file.close()
            os.remove(os.path.join(os.getcwd(),'texer.log'))
        else:
            print 'There were err…