PageRenderTime 51ms CodeModel.GetById 16ms RepoModel.GetById 0ms app.codeStats 0ms

/pymaclab/dsge/solvers/modsolvers.py

https://github.com/TomAugspurger/pymaclab
Python | 2924 lines | 2868 code | 3 blank | 53 comment | 15 complexity | 95c7bd14fdd96ba23a8aac13e0e13496 MD5 | raw file
Possible License(s): Apache-2.0

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

    

  1. '''

    2. 

  2. .. module:: modsolvers

    3. 

  3.    :platform: Linux

    4. 

  4.    :synopsis: The modsolvers module is an integral part of PyMacLab and is called extensively from the DSGEmodel class in module

    5. 

  5.               macrolab. It contains all of the dynamic solver methods most of which make use of the DSGE models' derivatives, i.e.

    6. 

  6.               the computed (numerical) Jacobian and Hessian.

    7. 

  7. 

    8. 

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

    9. 

  9. 

    10. 

  10. 

    11. 

  11. '''

    12. 

  12. 

    13. 

  13. from numpy import matlib as MAT

    14. 

  14. from numpy import linalg as LIN

    15. 

  15. from numpy.linalg import matrix_rank

    16. 

  16. from .. import helpers as HLP

    17. 

  17. from scipy.linalg import qz

    18. 

  18. import numpy as np

    19. 

  19. import pylab as P

    20. 

  20. from matplotlib import pyplot as PLT

    21. 

  21. import copy as COP

    22. 

  22. from pymaclab.filters import hpfilter

    23. 

  23. from pymaclab.filters import bkfilter

    24. 

  24. from pymaclab.filters import cffilter

    25. 

  25. from isolab import isolab

    26. 

  26. try:

    27. 

  27.     import mlabraw

    28. 

  28. except:

    29. 

  29.     pass

    30. 

  30. 

    31. 

  31. 

    32. 

  32. class MODsolvers(object):

    33. 

  33.     def __init__(self):

    34. 

  34.         pass

    35. 

  35. 

    36. 

  36. class Dynare(object):

    37. 

  37.     

    38. 

  38.     def __init__(self,attdic):

    39. 

  39.         self.attdic = attdic

    40. 

  40.         for keyo in attdic.keys():

    41. 

  41.             self.__setattr__(keyo,COP.deepcopy(attdic[keyo]))

    42. 

  42.     

    43. 

  43. 

    44. 

  44. class PyUhlig(MODsolvers):

    45. 

  45.     def __init__(self,intup):

    46. 

  46.         self.ntup = intup[0]

    47. 

  47.         self.nendo = self.ntup[0]

    48. 

  48.         self.ncon = self.ntup[1]

    49. 

  49.         self.nexo = self.ntup[2]

    50. 

  50. 

    51. 

  51.         self.eqindx = intup[1]

    52. 

  52.         self.vreg = intup[2]

    53. 

  53.         self.llsys_list = intup[3]

    54. 

  54.         self.diffli1 = intup[4]

    55. 

  55.         self.diffli2 = intup[5]

    56. 

  56.         diffli1 = self.diffli1

    57. 

  57.         diffli2 = self.diffli2

    58. 

  58. 

    59. 

  59.         deteq = []

    60. 

  60.         for x in self.eqindx['det']:

    61. 

  61.             deteq.append(self.llsys_list[x])

    62. 

  62.         self.deteq = deteq

    63. 

  63.         expeq = []

    64. 

  64.         for x in self.eqindx['exp']:

    65. 

  65.             expeq.append(self.llsys_list[x])

    66. 

  66.         self.expeq = expeq

    67. 

  67.         erreq = []

    68. 

  68.         for x in self.eqindx['err']:

    69. 

  69.             erreq.append(self.llsys_list[x])

    70. 

  70.         self.erreq = erreq

    71. 

  71. 

    72. 

  72.         detdif1 = []

    73. 

  73.         detdif2 = []

    74. 

  74.         for x in self.eqindx['det']:

    75. 

  75.             detdif1.append(diffli1[x])

    76. 

  76.             detdif2.append(diffli2[x])

    77. 

  77.         expdif1 = []

    78. 

  78.         expdif2 = []

    79. 

  79.         for x in self.eqindx['exp']:

    80. 

  80.             expdif1.append(diffli1[x])

    81. 

  81.             expdif2.append(diffli2[x])

    82. 

  82.         errdif1 = []

    83. 

  83.         errdif2 = []

    84. 

  84.         for x in self.eqindx['err']:

    85. 

  85.             errdif1.append(diffli1[x])

    86. 

  86.             errdif2.append(diffli2[x])

    87. 

  87.         self.detdif1 = detdif1

    88. 

  88.         self.detdif2 = detdif2

    89. 

  89.         self.expdif1 = expdif1

    90. 

  90.         self.expdif2 = expdif2

    91. 

  91.         self.errdif1 = errdif1

    92. 

  92.         self.errdif2 = errdif2

    93. 

  93. 

    94. 

  94.         self.mkmats()

    95. 

  95. 

    96. 

  96.     def mkmats(self):

    97. 

  97.         deteq = self.deteq

    98. 

  98.         expeq = self.expeq

    99. 

  99.         erreq = self.erreq

    100. 

  100.         nendo = self.nendo

    101. 

  101.         ncon = self.ncon

    102. 

  102.         nexo = self.nexo        

    103. 

  103.         vreg = self.vreg

    104. 

  104.         detdif1 = self.detdif1

    105. 

  105.         detdif2 = self.detdif2

    106. 

  106.         expdif1 = self.expdif1

    107. 

  107.         expdif2 = self.expdif2

    108. 

  108.         errdif1 = self.errdif1

    109. 

  109.         errdif2 = self.errdif2

    110. 

  110. 

    111. 

  111.         # Create argument list for matrix creation

    112. 

  112.         argli = [['AA',(detdif2,detdif1),deteq,(len(deteq),nendo),(None,'endo','0')],

    113. 

  113.              ['BB',(detdif2,detdif1),deteq,(len(deteq),nendo),(None,'endo','-1')],

    114. 

  114.              ['CC',(detdif2,detdif1),deteq,(len(deteq),ncon),(None,'con','0')],

    115. 

  115.              ['DD',(detdif2,detdif1),deteq,(len(deteq),nexo),(None,'exo','0')],

    116. 

  116.              ['FF',(expdif2,expdif1),expeq,(len(expeq),nendo),('0','endo','1')],

    117. 

  117.              ['GG',(expdif2,expdif1),expeq,(len(expeq),nendo),(None,'endo','0')],

    118. 

  118.              ['HH',(expdif2,expdif1),expeq,(len(expeq),nendo),(None,'endo','-1')],

    119. 

  119.              ['JJ',(expdif2,expdif1),expeq,(len(expeq),ncon),('0','con','1')],

    120. 

  120.              ['KK',(expdif2,expdif1),expeq,(len(expeq),ncon),(None,'con','0')],

    121. 

  121.              ['LL',(expdif2,expdif1),expeq,(len(expeq),nexo),('0','exo','1')],

    122. 

  122.              ['MM',(expdif2,expdif1),expeq,(len(expeq),nexo),(None,'exo','0')],

    123. 

  123.              ['NN',(errdif2,errdif1),erreq,(len(erreq),nexo),(None,'exo','0')]]

    124. 

  124. 

    125. 

  125.         # Create all matrices in a loop over argli

    126. 

  126.         for argx in argli:

    127. 

  127.             XX = MAT.zeros(argx[3])

    128. 

  128.             dic1 = dict([[x,'0'] for x in range(argx[3][1])])

    129. 

  129.             sXX = dict([[x,COP.deepcopy(dic1)] for x in range(len(argx[2]))])

    130. 

  130.             for y,x in enumerate(argx[2]):

    131. 

  131.                 if vreg(argx[4],x,False,'max'):

    132. 

  132.                     cont=[[z[0],z[1][1]] for z in vreg(argx[4],x,True,'min')]

    133. 

  133.                     for z in cont:

    134. 

  134.                         XX[y,z[1]] = argx[1][0][y][z[0]]

    135. 

  135.                         sXX[y][z[1]] = argx[1][1][y][z[0]]

    136. 

  136.             exec('self.'+argx[0]+' = XX')

    137. 

  137.             exec('self.'+'s'+argx[0]+' = sXX')

    138. 

  138. 

    139. 

  139.     def solve(self,sol_method='do_QZ',Xi_select='all'):

    140. 

  140.         self.output = {}

    141. 

  141.         self.Xi_select = Xi_select

    142. 

  142.         self.sol_method = sol_method

    143. 

  143.         if self.sol_method == 'do_PP':

    144. 

  144.             self.do_PP(self.Xi_select)

    145. 

  145.         elif self.sol_method == 'do_QZ':

    146. 

  146.             self.do_QZ(self.Xi_select,Tol=1.0000e-06)

    147. 

  147.         self.do_QRS()

    148. 

  148. 

    149. 

  149.     def do_QZ(self,Xi_select='all',Tol=1.0000e-06):

    150. 

  150.         # Make uhlig matrices locally available for computations

    151. 

  151.         AA = self.AA

    152. 

  152.         BB = self.BB

    153. 

  153.         CC = self.CC

    154. 

  154.         DD = self.DD

    155. 

  155.         FF = self.FF

    156. 

  156.         GG = self.GG

    157. 

  157.         HH = self.HH

    158. 

  158.         JJ = self.JJ

    159. 

  159.         KK = self.KK

    160. 

  160.         LL = self.LL

    161. 

  161.         MM = self.MM

    162. 

  162.         NN = self.NN

    163. 

  163.         q_expectational_equ = np.shape(FF)[0]

    164. 

  164.         m_states = np.shape(FF)[1]

    165. 

  165.         l_equ = np.shape(CC)[0]

    166. 

  166.         n_endog = np.shape(CC)[1]

    167. 

  167.         k_exog = min(np.shape(NN))

    168. 

  168. 

    169. 

  169. #        if HLP.rank(CC) < n_endog:

    170. 

  170.         if matrix_rank(CC) < n_endog:

    171. 

  171.             raise uhlerr, 'uhlerror: Rank(CC) needs to be at least\n\

    172. 

  172.                   equal to the number of endogenous variables'

    173. 

  173.         if l_equ == 0:

    174. 

  174.             CC_plus = LIN.pinv(CC)

    175. 

  175.             CC_0 = (HLP.null(CC.T)).T

    176. 

  176.             Psi_mat = FF

    177. 

  177.             Gamma_mat = -GG

    178. 

  178.             Theta_mat = -HH

    179. 

  179.             Xi_mat = \

    180. 

  180.                    np.concatenate((Gamma_mat,Theta_mat,MAT.eye(m_states),\

    181. 

  181.                           MAT.zeros((m_states,m_states))),1)

    182. 

  182.             Delta_mat = \

    183. 

  183.                   np.concatenate((Psi_mat,MAT.zeros((m_states,m_states)),\

    184. 

  184.                          MAT.zeros((m_states,m_states)),\

    185. 

  185.                          MAT.zeros((m_states,m_states)),MAT.eye(m_states),1))

    186. 

  186. 

    187. 

  187. 

    188. 

  188.         else:

    189. 

  189.             CC_plus = LIN.pinv(CC)

    190. 

  190.             CC_0 = HLP.null(CC.T)

    191. 

  191.             if l_equ > n_endog:

    192. 

  192.                 Psi_mat = \

    193. 

  193.                     np.concatenate((MAT.zeros((l_equ-n_endog,m_states)),FF-JJ*CC_plus*AA),1)

    194. 

  194.             elif l_equ == n_endog:

    195. 

  195.                 Psi_mat = np.concatenate((FF-JJ*CC_plus*AA),1)

    196. 

  196.             Gamma_mat = np.concatenate((CC_0*AA, JJ*CC_plus*BB-GG+KK*CC_plus*AA))

    197. 

  197.             Theta_mat = np.concatenate((CC_0*BB, KK*CC_plus*BB-HH))

    198. 

  198.             Xi_mat = \

    199. 

  199.                    np.concatenate((\

    200. 

  200.                        np.concatenate((Gamma_mat,Theta_mat),1),\

    201. 

  201.                        np.concatenate((MAT.identity(m_states),MAT.zeros((m_states,m_states))),1)\

    202. 

  202.                    ))

    203. 

  203.             Psi_mat = Psi_mat.reshape(m_states,m_states)

    204. 

  204.             Delta_mat = \

    205. 

  205.                   np.concatenate((\

    206. 

  206.                       np.concatenate((Psi_mat,MAT.zeros((m_states,m_states))),1),\

    207. 

  207.                       np.concatenate((MAT.zeros((m_states,m_states)),MAT.identity(m_states)),1)\

    208. 

  208.                   ))

    209. 

  209. 

    210. 

  210. #        (Delta_up,Xi_up,UUU,VVV)=\

    211. 

  211. #         HLP.qz(Delta_mat,Xi_mat,\

    212. 

  212. #            mode='complex',\

    213. 

  213. #            lapackname=lapackname,\

    214. 

  214. #            lapackpath=lapackpath)

    215. 

  215.         (Delta_up,Xi_up,UUU,VVV) = qz(Delta_mat,Xi_mat)

    216. 

  216.         

    217. 

  217. 

    218. 

  218.         d_Delta_up = MAT.diag(Delta_up)

    219. 

  219.         d_Xi_up = MAT.diag(Xi_up)

    220. 

  220.         Xi_eigval = MAT.zeros(N.shape(Delta_up))

    221. 

  221.         for i1 in range(0,N.shape(Delta_up)[0],1):

    222. 

  222.             Xi_eigval[i1,i1] = d_Xi_up[i1]/d_Delta_up[i1]

    223. 

  223.         d_Xi_eigval = np.diag(Xi_eigval)

    224. 

  224.         mat_tmp = MAT.zeros((N.shape(Xi_eigval)[0],3))

    225. 

  225.         i1=0

    226. 

  226.         for x in d_Xi_eigval:

    227. 

  227.             mat_tmp[i1,0] = d_Xi_eigval[i1]

    228. 

  228.             mat_tmp[i1,1] = abs(d_Xi_eigval)[i1]

    229. 

  229.             mat_tmp[i1,2] = i1

    230. 

  230.             i1=i1+1

    231. 

  231. 

    232. 

  232.         Xi_sortval = HLP.sortrows(mat_tmp,1)[:,0]

    233. 

  233.         # Need to do an argsort() on index column to turn to integers (Booleans) !!

    234. 

  234.         Xi_sortindex = HLP.sortrows(mat_tmp,1)[:,2].argsort(axis=0)

    235. 

  235.         Xi_sortabs = HLP.sortrows(mat_tmp,1)[:,1]

    236. 

  236. 

    237. 

  237.         # Root selection branch with Xi_select

    238. 

  238.         if Xi_select == 'all':

    239. 

  239.             Xi_select = np.arange(0,m_states,1)

    240. 

  240. 

    241. 

  241.         stake = max(abs(Xi_sortval[Xi_select])) + Tol

    242. 

  242.         (Delta_up,Xi_up,UUU,VVV) = self.qzdiv(stake,Delta_up,Xi_up,UUU,VVV)

    243. 

  243.         Lamda_mat = np.diag(Xi_sortval[Xi_select])

    244. 

  244.         trVVV = VVV.conjugate().T

    245. 

  245.         VVV_2_1 = trVVV[m_states:2*m_states,0:m_states]

    246. 

  246.         VVV_2_2 = trVVV[m_states:2*m_states,m_states:2*m_states]

    247. 

  247.         UUU_2_1 = UUU[m_states:2*m_states,0:m_states]

    248. 

  248. 

    249. 

  249.         PP = -(VVV_2_1.I*VVV_2_2)

    250. 

  250.         PP = np.real(PP)

    251. 

  251. 

    252. 

  252.         self.PP = PP

    253. 

  253.         self.CC_plus = CC_plus

    254. 

  254. 

    255. 

  255.     def do_PP(self,Xi_select='all'):

    256. 

  256.         # Make uhlig matrices locally available for computations

    257. 

  257.         AA = self.AA

    258. 

  258.         BB = self.BB

    259. 

  259.         CC = self.CC

    260. 

  260.         DD = self.DD

    261. 

  261.         FF = self.FF

    262. 

  262.         GG = self.GG

    263. 

  263.         HH = self.HH

    264. 

  264.         JJ = self.JJ

    265. 

  265.         KK = self.KK

    266. 

  266.         LL = self.LL

    267. 

  267.         MM = self.MM

    268. 

  268.         NN = self.NN

    269. 

  269.         q_expectational_equ = np.shape(FF)[0]

    270. 

  270.         m_states = np.shape(FF)[1]

    271. 

  271.         l_equ = np.shape(CC)[0]

    272. 

  272.         n_endog = np.shape(CC)[1]

    273. 

  273.         k_exog = min(N.shape(NN))

    274. 

  274. 

    275. 

  275. #        if HLP.rank(CC) < n_endog:

    276. 

  276.         if matrix_rank(CC) < n_endog:

    277. 

  277.             raise uhlerr, 'uhlerror: Rank(CC) needs to be at least\n\

    278. 

  278.                   equal to the number of endogenous variables'

    279. 

  279.         if l_equ == 0:

    280. 

  280.             CC_plus = LIN.pinv(CC)

    281. 

  281.             CC_0 = (HLP.null(CC.T)).T

    282. 

  282.             Psi_mat = FF

    283. 

  283.             Gamma_mat = -GG

    284. 

  284.             Theta_mat = -HH

    285. 

  285.             Xi_mat = \

    286. 

  286.                    np.concatenate((Gamma_mat,Theta_mat,MAT.eye(m_states),\

    287. 

  287.                           MAT.zeros((m_states,m_states))),1)

    288. 

  288.             Delta_mat = \

    289. 

  289.                   np.concatenate((Psi_mat,MAT.zeros((m_states,m_states)),\

    290. 

  290.                          MAT.zeros((m_states,m_states)),\

    291. 

  291.                          MAT.zeros((m_states,m_states)),MAT.eye(m_states),1))

    292. 

  292. 

    293. 

  293. 

    294. 

  294.         else:

    295. 

  295.             CC_plus = LIN.pinv(CC)

    296. 

  296.             CC_0 = HLP.null(CC.T)

    297. 

  297.             if l_equ > n_endog:

    298. 

  298.                 Psi_mat = \

    299. 

  299.                     np.concatenate((MAT.zeros((l_equ-n_endog,m_states)),FF-JJ*CC_plus*AA),1)

    300. 

  300.             elif l_equ == n_endog:

    301. 

  301.                 Psi_mat = np.concatenate((FF-JJ*CC_plus*AA),1)

    302. 

  302.             Gamma_mat = np.concatenate((CC_0*AA, JJ*CC_plus*BB-GG+KK*CC_plus*AA))

    303. 

  303.             Theta_mat = np.concatenate((CC_0*BB, KK*CC_plus*BB-HH))

    304. 

  304.             Xi_mat = \

    305. 

  305.                    np.concatenate((\

    306. 

  306.                        np.concatenate((Gamma_mat,Theta_mat),1),\

    307. 

  307.                        np.concatenate((MAT.eye(m_states),MAT.zeros((m_states,m_states))),1)\

    308. 

  308.                    ))

    309. 

  309.             Delta_mat = \

    310. 

  310.                   np.concatenate((\

    311. 

  311.                       np.concatenate((Psi_mat,MAT.zeros((m_states,m_states))),1),\

    312. 

  312.                       np.concatenate((MAT.zeros((m_states,m_states)),MAT.eye(m_states)),1)\

    313. 

  313.                   ))

    314. 

  314. 

    315. 

  315.         (Xi_eigvec,Xi_eigval) = HLP.eig(Xi_mat,Delta_mat)

    316. 

  316.         tmp_mat = MAT.eye(N.rank(Xi_eigvec))

    317. 

  317.         for i1 in range(0,N.rank(Xi_eigvec),1):

    318. 

  318.             tmp_mat[i1,i1] = float(Xi_eigval[0,i1])

    319. 

  319.         Xi_eigval = tmp_mat

    320. 

  320. 

    321. 

  321. 

    322. 

  322. #        if HLP.rank(Xi_eigvec) < m_states:

    323. 

  323.         if matrix_rank(Xi_eigvec) < m_states:

    324. 

  324.             raise uhlerr, 'uhlerror: Xi_mat is not diagonalizable!\n\

    325. 

  325.                   Cannot solve for PP. Maybe you should try the Schur-Decomposition\n\

    326. 

  326.                   method instead, use do_QZ()!!'

    327. 

  327. 

    328. 

  328. 

    329. 

  329.         d_Xi_eigval = np.diag(Xi_eigval)

    330. 

  330.         mat_tmp = MAT.zeros((N.rank(Xi_eigval),3))

    331. 

  331.         i1=0

    332. 

  332.         for x in d_Xi_eigval:

    333. 

  333.             mat_tmp[i1,0] = d_Xi_eigval[i1]

    334. 

  334.             mat_tmp[i1,1] = abs(d_Xi_eigval[i1])

    335. 

  335.             mat_tmp[i1,2] = i1

    336. 

  336.             i1=i1+1

    337. 

  337.         Xi_sortval = HLP.sortrows(mat_tmp,1)[:,0]

    338. 

  338.         # Need to do an argsort() on index column to turn to integers (Booleans) !!

    339. 

  339.         Xi_sortindex = HLP.sortrows(mat_tmp,1)[:,2].argsort(axis=0)

    340. 

  340.         Xi_sortabs = HLP.sortrows(mat_tmp,1)[:,1]

    341. 

  341.         Xi_sortvec = Xi_eigvec[0:2*m_states,:].take(Xi_sortindex.T.A, axis=1)

    342. 

  342. 

    343. 

  343.         # Root selection branch with Xi_select

    344. 

  344.         if Xi_select == 'all':

    345. 

  345.             Xi_select = np.arange(0,m_states,1)

    346. 

  346. 

    347. 

  347.         #Create Lambda_mat and Omega_mat

    348. 

  348.         Lambda_mat = MAT.zeros((len(Xi_select),len(Xi_select)))

    349. 

  349.         for i1 in range(0,len(Xi_select),1):

    350. 

  350.             Lambda_mat[i1,i1] = Xi_sortval[Xi_select][i1]

    351. 

  351.         Omega_mat = Xi_sortvec[m_states:2*m_states,:].take(Xi_select, axis=1)

    352. 

  352. 

    353. 

  353. #        if HLP.rank(Omega_mat) < m_states:

    354. 

  354.         if matrix_rank(Omega_mat) < m_states:

    355. 

  355.             raise uhlerr, 'uhlerror: Omega_mat is not invertible!!\n\

    356. 

  356.                   Therefore, solution for PP is not available'

    357. 

  357. 

    358. 

  358.         PP = Omega_mat*HLP.ctr((HLP.ctr(Omega_mat).I*HLP.ctr(Lambda_mat)))

    359. 

  359.         self.PP = PP

    360. 

  360.         self.CC_plus = CC_plus

    361. 

  361. 

    362. 

  362.     def do_QRS(self):

    363. 

  363.         # Make uhlig matrices locally available for computations

    364. 

  364.         AA = self.AA

    365. 

  365.         BB = self.BB

    366. 

  366.         CC = self.CC

    367. 

  367.         DD = self.DD

    368. 

  368.         FF = self.FF

    369. 

  369.         GG = self.GG

    370. 

  370.         HH = self.HH

    371. 

  371.         JJ = self.JJ

    372. 

  372.         KK = self.KK

    373. 

  373.         LL = self.LL

    374. 

  374.         MM = self.MM

    375. 

  375.         NN = self.NN

    376. 

  376.         PP = self.PP

    377. 

  377.         CC_plus = self.CC_plus

    378. 

  378.         (l_equ,m_states) = MAT.shape(AA)

    379. 

  379.         (l_equ,n_endog) = MAT.shape(CC)

    380. 

  380.         (l_equ,k_exog) = MAT.shape(DD)

    381. 

  381.         if l_equ == 0:

    382. 

  382.             RR = MAT.zeros((0,m_states))

    383. 

  383.             VV1 = MAT.kron(MAT.transpose(NN),FF)+MAT.kron(MAT.identity(k_exog),FF*PP+GG)

    384. 

  384.             VV2 = MAT.kron(MAT.transpose(NN),JJ)+MAT.kron(MAT.identity(k_exog),KK)

    385. 

  385.             VV = MAT.hstack((VV1,VV2))

    386. 

  386.         else:

    387. 

  387.             RR = -CC_plus*(AA*PP+BB)

    388. 

  388.             VV1 = MAT.kron(MAT.identity(k_exog),AA)

    389. 

  389.             VV2 = MAT.kron(MAT.identity(k_exog),CC)

    390. 

  390.             VV3 = MAT.kron(MAT.transpose(NN),FF)+MAT.kron(MAT.identity(k_exog),FF*PP+JJ*RR+GG)

    391. 

  391.             VV4 = MAT.kron(MAT.transpose(NN),JJ)+MAT.kron(MAT.identity(k_exog),KK)

    392. 

  392.             VV = MAT.vstack((MAT.hstack((VV1,VV2)),MAT.hstack((VV3,VV4))))

    393. 

  393.         self.RR = RR

    394. 

  394. 

    395. 

  395.         LLNN_plus_MM = LL*NN+MM

    396. 

  396.         NON = MAT.hstack((DD.flatten(1),LLNN_plus_MM.flatten(1)))

    397. 

  397.         try:

    398. 

  398.             QQSS_vec = -(VV.I*MAT.transpose(NON))

    399. 

  399.         except MyErr:

    400. 

  400.             print 'Error: Matrix VV is not invertible!'

    401. 

  401.         QQ = QQSS_vec[0:m_states*k_exog,:].reshape(-1,m_states).transpose()

    402. 

  402.         SS = QQSS_vec[(m_states*k_exog):((m_states+n_endog)*k_exog),:].reshape(-1,n_endog).transpose()

    403. 

  403.         WW1 = MAT.hstack((MAT.identity(m_states),MAT.zeros((m_states,k_exog))))

    404. 

  404.         WW2 = MAT.hstack((RR*LIN.pinv(PP),SS-RR*LIN.pinv(PP)*QQ))

    405. 

  405.         WW3 = MAT.hstack((MAT.zeros((k_exog,m_states)),MAT.identity(k_exog)))

    406. 

  406.         WW = MAT.vstack((WW1,WW2,WW3))

    407. 

  407.         self.WW = WW

    408. 

  408.         self.QQ = QQ

    409. 

  409.         self.SS = SS

    410. 

  410.         del self.CC_plus

    411. 

  411. 

    412. 

  412.     def qzdiv(self,stake,A,B,Q,Z):

    413. 

  413.         n = np.shape(A)[0]

    414. 

  414.         root = np.hstack((abs(N.mat(N.diag(A)).T),abs(N.mat(N.diag(B)).T)))

    415. 

  415.         index_mat = (root[:,0]<1e-13).choose(root[:,0],1)

    416. 

  416.         index_mat = (index_mat>1e-13).choose(root[:,0],0)

    417. 

  417.         root[:,0] = root[:,0]-MAT.multiply(index_mat,(root[:,0]+root[:,1]))

    418. 

  418.         root[:,1] = root[:,1]/root[:,0]

    419. 

  419.         for i1 in range(n-1,-1,-1):

    420. 

  420.             m='none'

    421. 

  421.             for i2 in range(i1,-1,-1):

    422. 

  422.                 if root[i2,1] > stake or root[i2,1] < -0.1:

    423. 

  423.                     m=i2

    424. 

  424.                     break

    425. 

  425.             if m == 'none':

    426. 

  426.                 return (A,B,Q,Z)

    427. 

  427.             else:

    428. 

  428.                 for i3 in range(m,i1,1):

    429. 

  429.                     (A,B,Q,Z) = self.qzswitch(i3,A,B,Q,Z)

    430. 

  430.                     tmp = COP.deepcopy(root[i3,1])

    431. 

  431.                     root[i3,1] = root[i3+1,1]

    432. 

  432.                     root[i3+1,1] = tmp

    433. 

  433. 

    434. 

  434.     def qzswitch(self,i,A,B,Q,Z):

    435. 

  435.         a = A[i,i]

    436. 

  436.         d = B[i,i]

    437. 

  437.         b = A[i,i+1]

    438. 

  438.         e = B[i,i+1]

    439. 

  439.         c = A[i+1,i+1]

    440. 

  440.         f = B[i+1,i+1]

    441. 

  441.         wz = np.mat(N.hstack(((c*e-f*b),(c*d-f*a).conjugate().T)))

    442. 

  442.         xy = np.mat(N.hstack(((b*d-e*a).conjugate().T,(c*d-f*a).conjugate().T)))

    443. 

  443.         n = np.mat(N.sqrt(wz*wz.conjugate().T))

    444. 

  444.         m = np.mat(N.sqrt(xy*xy.conjugate().T))

    445. 

  445.         if n.all() == 0:

    446. 

  446.             return

    447. 

  447.         else:

    448. 

  448.             wz = np.mat(wz/n)

    449. 

  449.             xy = np.mat(xy/m)

    450. 

  450.             wz = np.vstack((wz,N.hstack((-wz[:,1].T,wz[:,0].T))))

    451. 

  451.             xy = np.vstack((xy,N.hstack((-xy[:,1].T,xy[:,0].T))))

    452. 

  452.             A[i:i+2,:] = xy*A[i:i+2,:]

    453. 

  453.             B[i:i+2,:] = xy*B[i:i+2,:]

    454. 

  454.             A[:,i:i+2] = A[:,i:i+2]*wz

    455. 

  455.             B[:,i:i+2] = B[:,i:i+2]*wz

    456. 

  456.             Z[:,i:i+2] = Z[:,i:i+2]*wz

    457. 

  457.             Q[i:i+2,:] = xy*Q[i:i+2,:]

    458. 

  458.         return (A,B,Q,Z)

    459. 

  459. 

    460. 

  460. #----------------------------------------------------------------------------------------------------------------------

    461. 

  461. class MatUhlig:

    462. 

  462. 

    463. 

  463.     def __init__(self,intup):

    464. 

  464.         self.ntup = intup[0]

    465. 

  465.         self.nendo = self.ntup[0]

    466. 

  466.         self.ncon = self.ntup[1]

    467. 

  467.         self.nexo = self.ntup[2]

    468. 

  468. 

    469. 

  469.         self.eqindx = intup[1]

    470. 

  470.         self.vreg = intup[2]

    471. 

  471.         self.llsys_list = intup[3]

    472. 

  472.         self.diffli1 = intup[4]

    473. 

  473.         self.diffli2 = intup[5]

    474. 

  474.         self.sess1 = intup[6]

    475. 

  475.         self.vardic = intup[7]

    476. 

  476.         diffli1 = self.diffli1

    477. 

  477.         diffli2 = self.diffli2

    478. 

  478. 

    479. 

  479.         deteq = []

    480. 

  480.         for x in self.eqindx['det']:

    481. 

  481.             deteq.append(self.llsys_list[x])

    482. 

  482.         self.deteq = deteq

    483. 

  483.         expeq = []

    484. 

  484.         for x in self.eqindx['exp']:

    485. 

  485.             expeq.append(self.llsys_list[x])

    486. 

  486.         self.expeq = expeq

    487. 

  487.         erreq = []

    488. 

  488.         for x in self.eqindx['err']:

    489. 

  489.             erreq.append(self.llsys_list[x])

    490. 

  490.         self.erreq = erreq

    491. 

  491. 

    492. 

  492.         detdif1 = []

    493. 

  493.         detdif2 = []

    494. 

  494.         for x in self.eqindx['det']:

    495. 

  495.             detdif1.append(diffli1[x])

    496. 

  496.             detdif2.append(diffli2[x])

    497. 

  497.         expdif1 = []

    498. 

  498.         expdif2 = []

    499. 

  499.         for x in self.eqindx['exp']:

    500. 

  500.             expdif1.append(diffli1[x])

    501. 

  501.             expdif2.append(diffli2[x])

    502. 

  502.         errdif1 = []

    503. 

  503.         errdif2 = []

    504. 

  504.         for x in self.eqindx['err']:

    505. 

  505.             errdif1.append(diffli1[x])

    506. 

  506.             errdif2.append(diffli2[x])

    507. 

  507.         self.detdif1 = detdif1

    508. 

  508.         self.detdif2 = detdif2

    509. 

  509.         self.expdif1 = expdif1

    510. 

  510.         self.expdif2 = expdif2

    511. 

  511.         self.errdif1 = errdif1

    512. 

  512.         self.errdif2 = errdif2

    513. 

  513. 

    514. 

  514.         self.mkmats()

    515. 

  515. 

    516. 

  516.     def mkmats(self):

    517. 

  517.         deteq = self.deteq

    518. 

  518.         expeq = self.expeq

    519. 

  519.         erreq = self.erreq

    520. 

  520.         nendo = self.nendo

    521. 

  521.         ncon = self.ncon

    522. 

  522.         nexo = self.nexo        

    523. 

  523.         vreg = self.vreg

    524. 

  524.         detdif1 = self.detdif1

    525. 

  525.         detdif2 = self.detdif2

    526. 

  526.         expdif1 = self.expdif1

    527. 

  527.         expdif2 = self.expdif2

    528. 

  528.         errdif1 = self.errdif1

    529. 

  529.         errdif2 = self.errdif2

    530. 

  530. 

    531. 

  531.         # Create argument list for matrix creation

    532. 

  532.         argli = [['AA',(detdif2,detdif1),deteq,(len(deteq),nendo),(None,'endo','0')],

    533. 

  533.              ['BB',(detdif2,detdif1),deteq,(len(deteq),nendo),(None,'endo','-1')],

    534. 

  534.              ['CC',(detdif2,detdif1),deteq,(len(deteq),ncon),(None,'con','0')],

    535. 

  535.              ['DD',(detdif2,detdif1),deteq,(len(deteq),nexo),(None,'exo','0')],

    536. 

  536.              ['FF',(expdif2,expdif1),expeq,(len(expeq),nendo),('0','endo','1')],

    537. 

  537.              ['GG',(expdif2,expdif1),expeq,(len(expeq),nendo),(None,'endo','0')],

    538. 

  538.              ['HH',(expdif2,expdif1),expeq,(len(expeq),nendo),(None,'endo','-1')],

    539. 

  539.              ['JJ',(expdif2,expdif1),expeq,(len(expeq),ncon),('0','con','1')],

    540. 

  540.              ['KK',(expdif2,expdif1),expeq,(len(expeq),ncon),(None,'con','0')],

    541. 

  541.              ['LL',(expdif2,expdif1),expeq,(len(expeq),nexo),('0','exo','1')],

    542. 

  542.              ['MM',(expdif2,expdif1),expeq,(len(expeq),nexo),(None,'exo','0')],

    543. 

  543.              ['NN',(errdif2,errdif1),erreq,(len(erreq),nexo),(None,'exo','0')]]

    544. 

  544. 

    545. 

  545.         # Create all matrices in a loop over argli

    546. 

  546.         for argx in argli:

    547. 

  547.             XX = MAT.zeros(argx[3])

    548. 

  548.             dic1 = dict([[x,'0'] for x in range(argx[3][1])])

    549. 

  549.             sXX = dict([[x,COP.deepcopy(dic1)] for x in range(len(argx[2]))])

    550. 

  550.             for y,x in enumerate(argx[2]):

    551. 

  551.                 if vreg(argx[4],x,False,'max'):

    552. 

  552.                     cont=[[z[0],z[1][1]] for z in vreg(argx[4],x,True,'min')]

    553. 

  553.                     for z in cont:

    554. 

  554.                         XX[y,z[1]] = argx[1][0][y][z[0]]

    555. 

  555.                         sXX[y][z[1]] = argx[1][1][y][z[0]]

    556. 

  556.             exec('self.'+argx[0]+' = XX')

    557. 

  557.             exec('self.'+'s'+argx[0]+' = sXX')

    558. 

  558. 

    559. 

  559.     def solve(self,sol_method='do_QZ',Xi_select='all'):

    560. 

  560.         # Create varnames with correct string length

    561. 

  561.         tmp_list =[x[1] for x in self.vardic['endo']['var']]\

    562. 

  562.              +[x[1] for x in self.vardic['con']['var']]\

    563. 

  563.              +[x[1] for x in self.vardic['exo']['var']]

    564. 

  564.         tmp_list2 = [[len(x),x] for x in tmp_list]

    565. 

  565.         tmp_list2.sort()

    566. 

  566.         max_len = tmp_list2[-1][0]

    567. 

  567.         i1=0

    568. 

  568.         for x in tmp_list:

    569. 

  569.             tmp_list[i1]=x+(max_len-len(x))*' '

    570. 

  570.             i1=i1+1

    571. 

  571.         varnstring = 'VARNAMES = ['

    572. 

  572.         for x in tmp_list:

    573. 

  573.             varnstring = varnstring + "'" +x[1]+ "'" + ','

    574. 

  574.         varnstring = varnstring[0:-1]+'];'

    575. 

  575. 

    576. 

  576.         # Start matlab session and calculate

    577. 

  577.         sess1 = self.sess1

    578. 

  578.         mlabraw.eval(sess1,'clear all;')

    579. 

  579.         mlabraw.eval(sess1,varnstring)

    580. 

  580.         mlabraw.put(sess1,'AA',self.AA)

    581. 

  581.         mlabraw.put(sess1,'BB',self.BB)

    582. 

  582.         mlabraw.put(sess1,'CC',self.CC)

    583. 

  583.         mlabraw.put(sess1,'DD',self.DD)

    584. 

  584.         mlabraw.put(sess1,'FF',self.FF)

    585. 

  585.         mlabraw.put(sess1,'GG',self.GG)

    586. 

  586.         mlabraw.put(sess1,'HH',self.HH)

    587. 

  587.         mlabraw.put(sess1,'JJ',self.JJ)

    588. 

  588.         mlabraw.put(sess1,'KK',self.KK)

    589. 

  589.         mlabraw.put(sess1,'LL',self.LL)

    590. 

  590.         mlabraw.put(sess1,'MM',self.MM)

    591. 

  591.         mlabraw.put(sess1,'NN',self.NN)

    592. 

  592.         mlabraw.eval(sess1,'cd '+mlabpath)

    593. 

  593.         mlabraw.eval(sess1,'cd Toolkit41')

    594. 

  594.         message = ' '*70

    595. 

  595.         eval_list = ['message='+"'"+message+"'",

    596. 

  596.                  'warnings = [];',

    597. 

  597.                  'options;',

    598. 

  598.                  'solve;']

    599. 

  599.         try:

    600. 

  600.             for x in eval_list:

    601. 

  601.                 mlabraw.eval(sess1,x)

    602. 

  602.             self.PP = np.matrix(mlabraw.get(sess1,'PP'))

    603. 

  603.             self.QQ = np.matrix(mlabraw.get(sess1,'QQ'))

    604. 

  604.             self.RR = np.matrix(mlabraw.get(sess1,'RR'))

    605. 

  605.             self.SS = np.matrix(mlabraw.get(sess1,'SS'))

    606. 

  606.             self.WW = np.matrix(mlabraw.get(sess1,'WW'))

    607. 

  607.         except MyErr:

    608. 

  608.             pass

    609. 

  609.         finally:

    610. 

  610.             return

    611. 

  611. #----------------------------------------------------------------------------------------------------------------------

    612. 

  612. class MatKlein:

    613. 

  613. 

    614. 

  614.     def __init__(self,intup):

    615. 

  615.         self.sess1 = intup[-1]

    616. 

  616.         self.uhlig = PyUhlig(intup[:-1])

    617. 

  617.         self.uhlig.mkmats()

    618. 

  618.         self.AA = self.uhlig.AA

    619. 

  619.         self.BB = self.uhlig.BB

    620. 

  620.         self.CC = self.uhlig.CC

    621. 

  621.         self.DD = self.uhlig.DD

    622. 

  622.         self.FF = self.uhlig.FF

    623. 

  623.         self.GG = self.uhlig.GG

    624. 

  624.         self.HH = self.uhlig.HH

    625. 

  625.         self.JJ = self.uhlig.JJ

    626. 

  626.         self.KK = self.uhlig.KK

    627. 

  627.         self.LL = self.uhlig.LL

    628. 

  628.         self.MM = self.uhlig.MM

    629. 

  629.         self.NN = self.uhlig.NN

    630. 

  630.         self.mkKleinMats()

    631. 

  631.         self.outdic = {}

    632. 

  632. 

    633. 

  633.     def mkKleinMats(self):

    634. 

  634.         # Make uhlig matrices locally available for computations

    635. 

  635.         AA = self.AA

    636. 

  636.         BB = self.BB

    637. 

  637.         CC = self.CC

    638. 

  638.         DD = self.DD

    639. 

  639.         FF = self.FF

    640. 

  640.         GG = self.GG

    641. 

  641.         HH = self.HH

    642. 

  642.         JJ = self.JJ

    643. 

  643.         KK = self.KK

    644. 

  644.         LL = self.LL

    645. 

  645.         MM = self.MM

    646. 

  646.         NN = self.NN

    647. 

  647. 

    648. 

  648.         # Determine size of states, endogenous

    649. 

  649.         exo_st = MAT.shape(NN)[1]

    650. 

  650.         endo_st = MAT.shape(BB)[1]

    651. 

  651.         endo_cn = MAT.shape(CC)[1]

    652. 

  652.         n_deteq = MAT.shape(AA)[0]

    653. 

  653.         n_expeq = MAT.shape(JJ)[0]

    654. 

  654.         tot_st = exo_st+endo_st

    655. 

  655.         self.tstates = tot_st

    656. 

  656.         tot_var = tot_st+endo_cn

    657. 

  657. 

    658. 

  658.         klein_A_rtwo = MAT.hstack((LL,GG))

    659. 

  659.         klein_A_rtwo = MAT.hstack((klein_A_rtwo,JJ))

    660. 

  660.         klein_A_rtwo_rows = MAT.shape(klein_A_rtwo)[0]

    661. 

  661.         klein_A_rtwo_cols = MAT.shape(klein_A_rtwo)[1]

    662. 

  662. 

    663. 

  663.         klein_A_rone = MAT.zeros((exo_st,klein_A_rtwo_cols))

    664. 

  664.         klein_A_rone = MAT.hstack((MAT.identity(exo_st),klein_A_rone[:,exo_st:]))

    665. 

  665. 

    666. 

  666.         klein_A_rthree = MAT.hstack((MAT.zeros((n_deteq,exo_st)),AA))

    667. 

  667.         klein_A_rthree = MAT.hstack((klein_A_rthree,MAT.zeros((n_deteq,endo_cn))))

    668. 

  668. 

    669. 

  669.         klein_A = MAT.vstack((klein_A_rone,klein_A_rtwo))

    670. 

  670.         klein_A = MAT.vstack((klein_A,klein_A_rthree))

    671. 

  671. 

    672. 

  672.         klein_B_rone = MAT.zeros((exo_st,klein_A_rtwo_cols))

    673. 

  673. 

    674. 

  674.         klein_B_rone = MAT.hstack((NN,klein_B_rone[:,exo_st:]))

    675. 

  675.         klein_B_rtwo = MAT.hstack((-MM,-HH))

    676. 

  676.         klein_B_rtwo = MAT.hstack((klein_B_rtwo,-KK))

    677. 

  677.         klein_B_rthree = MAT.hstack((-DD,-BB))

    678. 

  678.         klein_B_rthree = MAT.hstack((klein_B_rthree,-CC))

    679. 

  679. 

    680. 

  680.         klein_B = MAT.vstack((klein_B_rone,klein_B_rtwo))

    681. 

  681.         klein_B = MAT.vstack((klein_B,klein_B_rthree))  

    682. 

  682. 

    683. 

  683.         self.A = klein_A

    684. 

  684.         self.B = klein_B

    685. 

  685. 

    686. 

  686.     def solve(self):

    687. 

  687.         A = self.A

    688. 

  688.         B = self.B

    689. 

  689.         tstates = self.tstates

    690. 

  690.         sess1 = self.sess1

    691. 

  691.         mlabraw.eval(sess1,'clear all;')

    692. 

  692.         mlabraw.eval(sess1,'cd '+mlabpath)

    693. 

  693.         mlabraw.eval(sess1,'cd Klein')

    694. 

  694.         mlabraw.put(sess1,'AA',A)

    695. 

  695.         mlabraw.put(sess1,'BB',B)

    696. 

  696.         mlabraw.put(sess1,'tstates',tstates)

    697. 

  697.         try:

    698. 

  698.             mlabraw.eval(sess1,'[F,P,Z11]=solab(AA,BB,tstates)')

    699. 

  699.             self.F = np.matrix(mlabraw.get(sess1,'F'))

    700. 

  700.             self.P = np.matrix(mlabraw.get(sess1,'P'))

    701. 

  701.             self.Z11 = np.matrix(mlabraw.get(sess1,'Z11'))

    702. 

  702.             self.outdic['P'] = self.P

    703. 

  703.             self.outdic['F'] = self.F

    704. 

  704.         except MyErr:

    705. 

  705.             pass

    706. 

  706.         finally:

    707. 

  707.             return

    708. 

  708. #----------------------------------------------------------------------------------------------------------------------

    709. 

  709. class MatKleinD:

    710. 

  710. 

    711. 

  711.     def __init__(self,intup):

    712. 

  712.         self.interm3 = intup[0]

    713. 

  713.         self.param = intup[1]

    714. 

  714.         self.sstate_list = intup[2]

    715. 

  715.         self.varvecs = intup[3]

    716. 

  716.         self.vardic = intup[4]

    717. 

  717.         self.vardic2 = intup[5]

    718. 

  718.         self.shockdic = intup[6]

    719. 

  719.         self.shockdic2 = intup[7]

    720. 

  720.         self.varns = intup[8]

    721. 

  721.         self.re_var = intup[9]

    722. 

  722.         self.re_var2 = intup[10]

    723. 

  723.         self.ishockdic = intup[11]

    724. 

  724.         self.ishockdic2 = intup[12]

    725. 

  725.         self.sess1 = intup[13]

    726. 

  726.         self.sstate = intup[14]

    727. 

  727.         self.mkeqs()

    728. 

  728.         self.mkvarl()

    729. 

  729.         self.mkeqs2()

    730. 

  730.         self.mkgrad()

    731. 

  731.         self.mkAB()

    732. 

  732.         self.solab()

    733. 

  733.         self.mkhess()

    734. 

  734.         self.solab2()

    735. 

  735. 

    736. 

  736.     def mkeqs(self):

    737. 

  737.         list_tmp = COP.deepcopy(self.interm3)

    738. 

  738.         reg2 = RE.compile('\*{2,2}')

    739. 

  739.         reva2 = self.re_var2

    740. 

  740.         i1=0

    741. 

  741.         for x in list_tmp:

    742. 

  742.             list_tmp[i1]=reg2.sub('^',x.split('=')[0]).strip()

    743. 

  743.             for y in self.subvars(tuple(self.varvecs['e'])):

    744. 

  744.                 while y in list_tmp[i1]:

    745. 

  745.                     list_tmp[i1] = list_tmp[i1].replace(y,'1')

    746. 

  746.             i1=i1+1

    747. 

  747.         self.interm4 = list_tmp

    748. 

  748. 

    749. 

  749.     def mkvarl(self):

    750. 

  750.         varlist=[]

    751. 

  751.         nexost = len(self.varvecs['z'])

    752. 

  752.         nendost = len(self.varvecs['k'])

    753. 

  753.         nendocon = len(self.varvecs['y'])

    754. 

  754.         exost_f = makeForward(self.varvecs['z'],'1','0')

    755. 

  755.         for x in exost_f:

    756. 

  756.             varlist.append(x)

    757. 

  757.         for x in self.varvecs['k']:

    758. 

  758.             varlist.append(x)

    759. 

  759.         endocons_f = makeForward(self.varvecs['y'],'1','0')

    760. 

  760.         for x in endocons_f:

    761. 

  761.             varlist.append(x)

    762. 

  762.         exost_l = self.varvecs['z']

    763. 

  763.         for x in exost_l:

    764. 

  764.             varlist.append(x)

    765. 

  765.         endost_l = makeLags(self.varvecs['k'],'1')

    766. 

  766.         for x in endost_l:

    767. 

  767.             varlist.append(x)

    768. 

  768.         for x in self.varvecs['y']:

    769. 

  769.             varlist.append(x)

    770. 

  770.         varlist2 = list(self.subvars(tuple(varlist)))

    771. 

  771.         self.vlist = varlist2

    772. 

  772. 

    773. 

  773.         varlist=[]

    774. 

  774.         for x in self.varvecs['z']:

    775. 

  775.             varlist.append(x)

    776. 

  776.         for x in self.varvecs['k']:

    777. 

  777.             varlist.append(x)

    778. 

  778.         endocons_f = makeForward(self.varvecs['y'],'1','0')

    779. 

  779.         for x in endocons_f:

    780. 

  780.             varlist.append(x)

    781. 

  781.         exost_l = makeLags(self.varvecs['z'],'1')

    782. 

  782.         for x in exost_l:

    783. 

  783.             varlist.append(x)

    784. 

  784.         endost_l = makeLags(self.varvecs['k'],'1')

    785. 

  785.         for x in endost_l:

    786. 

  786.             varlist.append(x)

    787. 

  787.         for x in self.varvecs['y']:

    788. 

  788.             varlist.append(x)

    789. 

  789.         varlist2 = list(self.subvars(tuple(varlist)))

    790. 

  790.         self.vlist2 = varlist2

    791. 

  791. 

    792. 

  792.     def mkeqs2(self):

    793. 

  793.         nexost = len(self.varvecs['z'])

    794. 

  794.         nendost = len(self.varvecs['k'])

    795. 

  795.         nendocon = len(self.varvecs['y'])

    796. 

  796.         list_tmp1 = COP.deepcopy(self.interm4)

    797. 

  797.         eq_len = len(self.interm4)

    798. 

  798.         sub_list2 = self.vlist2

    799. 

  799.         sub_list = self.vlist

    800. 

  800.         reva2 = self.re_var2

    801. 

  801.         i1=0

    802. 

  802.         for x in list_tmp1:

    803. 

  803.             str_tmp1 = x[:]

    804. 

  804.             i2=0

    805. 

  805.             for y1,y2 in zip(sub_list,sub_list2):

    806. 

  806.                 if i1 < eq_len-nexost:

    807. 

  807.                     while reva2.search(str_tmp1) != None:

    808. 

  808.                         ma1 = reva2.search(str_tmp1)

    809. 

  809.                         indx1 = sub_list.index(ma1.group(0))

    810. 

  810.                         str_tmp1 = str_tmp1.replace(ma1.group(0),'exp(x0('+str(indx1+1)+'))')[:]

    811. 

  811.                     i2 = i2 + 1

    812. 

  812.                 else:

    813. 

  813.                     while reva2.search(str_tmp1) != None:

    814. 

  814.                         ma2 = reva2.search(str_tmp1)

    815. 

  815.                         indx2 = sub_list2.index(ma2.group(0))

    816. 

  816.                         str_tmp1 = str_tmp1.replace(ma2.group(0),'exp(x0('+str(indx2+1)+'))')[:]

    817. 

  817.                     i2 = i2 + 1

    818. 

  818.             list_tmp1[i1] = str_tmp1

    819. 

  819.             i1=i1+1

    820. 

  820. 

    821. 

  821.         self.interm5 = list_tmp1        

    822. 

  822. 

    823. 

  823.     def mkfunc(self):

    824. 

  824.         rege = RE.compile('\*{2,2}')

    825. 

  825.         inde = ' '*2

    826. 

  826.         eq_len=len(self.interm5)

    827. 

  827.         mfile=open(path_opts['MlabPath']['path']+'/Klein/'+'mfunc.m','w')

    828. 

  828.         mfile.write('function y = mfunc(x0)\n')

    829. 

  829.         mfile.write(inde+'%Parameter Values\n')

    830. 

  830.         for x in self.param.items():

    831. 

  831.             mfile.write(inde+x[0]+'='+str(x[1])+';\n')

    832. 

  832.         mfile.write('\n')

    833. 

  833.         mfile.write(inde+'%Steady State Values\n')

    834. 

  834.         for x in self.sstate_list:

    835. 

  835.             x[0]=rege.sub('^',x[0])

    836. 

  836.             x[1]=rege.sub('^',x[1])

    837. 

  837.             mfile.write(inde+x[0]+'='+str(x[1])+';\n')

    838. 

  838.         mfile.write('\n')

    839. 

  839.         mfile.write(inde+'y=zeros('+str(eq_len)+',1);\n')

    840. 

  840.         mfile.write('\n')

    841. 

  841.         for i1,i2 in zip(range(1,eq_len+1,1),self.interm5):

    842. 

  842.             mfile.write(inde+'y('+str(i1)+')='+i2+';'+'\n')

    843. 

  843.         mfile.close()

    844. 

  844. 

    845. 

  845.     def mkgrad(self):

    846. 

  846.         self.mkfunc()

    847. 

  847.         sess1 = self.sess1

    848. 

  848.         sstate = self.sstate

    849. 

  849.         ssdic = sstate[0]

    850. 

  850.         ssdic.update(sstate[1])

    851. 

  851.         ssdic2 = {}

    852. 

  852.         for x in ssdic.keys():

    853. 

  853.             if x[-4:] == '_bar':

    854. 

  854.                 ssdic2[x] = ssdic[x]

    855. 

  855.         for x in ssdic2.keys():

    856. 

  856.             ssdic2[x] = np.log(ssdic2[x])

    857. 

  857.         vlist2 = []

    858. 

  858.         for x in self.vlist:

    859. 

  859.             ma = self.re_var2.search(x)

    860. 

  860.             tvar = ma.group('var')

    861. 

  861.             tvar = tvar.upper()

    862. 

  862.             tvar = tvar+'_bar'

    863. 

  863.             vlist2.append(tvar)

    864. 

  864.         invec = []

    865. 

  865.         for x in vlist2:

    866. 

  866.             invec.append(ssdic2[x])

    867. 

  867.         mlabraw.eval(sess1,'clear all;')

    868. 

  868.         mlabraw.eval(sess1,'cd '+path_opts['MlabPath']['path'])

    869. 

  869.         mlabraw.eval(sess1,'cd Klein')

    870. 

  870.         mlabraw.put(sess1,'x0',invec)

    871. 

  871.         mlabraw.eval(sess1,"grdm=centgrad('mfunc',x0)")

    872. 

  872.         gradm = np.matrix(mlabraw.get(sess1,'grdm'))

    873. 

  873.         self.gradm = gradm

    874. 

  874. 

    875. 

  875.     def mkAB(self):

    876. 

  876.         nexost = len(self.varvecs['z'])

    877. 

  877.         nendost = len(self.varvecs['k'])

    878. 

  878.         nendocon = len(self.varvecs['y'])

    879. 

  879.         A = self.gradm[:,:nexost+nendost+nendocon]

    880. 

  880.         B = -self.gradm[:,nexost+nendost+nendocon:]

    881. 

  881.         self.AA = A

    882. 

  882.         self.BB = B

    883. 

  883. 

    884. 

  884.     def solab(self):

    885. 

  885.         A = self.AA

    886. 

  886.         B = self.BB

    887. 

  887.         tstates = len(self.varvecs['k'])+len(self.varvecs['z'])

    888. 

  888.         (F,P,retcon) = isolab(A,B,tstates,MAT.shape(A)[0])

    889. 

  889.         if MAT.sum(P.reshape(-1,1)) == 0.0:

    890. 

  890.             return

    891. 

  891.         else:

    892. 

  892.             self.P = np.matrix(P)

    893. 

  893.             self.F = np.matrix(F)

    894. 

  894. 

    895. 

  895.     def mkhess(self):

    896. 

  896.         self.mkfunc()

    897. 

  897.         sess1 = self.sess1

    898. 

  898.         sstate = self.sstate

    899. 

  899.         ssdic = sstate[0]

    900. 

  900.         ssdic.update(sstate[1])

    901. 

  901.         ssdic2 = {}

    902. 

  902.         for x in ssdic.keys():

    903. 

  903.             if x[-4:] == '_bar':

    904. 

  904.                 ssdic2[x] = ssdic[x]

    905. 

  905.         for x in ssdic2.keys():

    906. 

  906.             ssdic2[x] = np.log(ssdic2[x])

    907. 

  907.         vlist2 = []

    908. 

  908.         for x in self.vlist:

    909. 

  909.             ma = self.re_var2.search(x)

    910. 

  910.             tvar = ma.group('var')

    911. 

  911.             tvar = tvar.upper()

    912. 

  912.             tvar = tvar+'_bar'

    913. 

  913.             vlist2.append(tvar)

    914. 

  914.         invec = []

    915. 

  915.         for x in vlist2:

    916. 

  916.             invec.append(ssdic2[x])

    917. 

  917.         mlabraw.eval(sess1,'clear all;')

    918. 

  918.         mlabraw.eval(sess1,'cd '+path_opts['MlabPath']['path'])

    919. 

  919.         mlabraw.eval(sess1,'cd Klein')

    920. 

  920.         mlabraw.put(sess1,'x0',invec)

    921. 

  921.         mlabraw.eval(sess1,"hessm=centhess('mfunc',x0)")

    922. 

  922.         hessm = np.matrix(mlabraw.get(sess1,'hessm'))

    923. 

  923.         self.hessm = hessm

    924. 

  924. 

    925. 

  925.     def solab2(self):

    926. 

  926.         sess1 = self.sess1

    927. 

  927.         mlabraw.eval(sess1,'clear all;')

    928. 

  928.         mlabraw.eval(sess1,'cd '+path_opts['MlabPath']['path'])

    929. 

  929.         mlabraw.eval(sess1,'cd Klein')

    930. 

  930.         mlabraw.put(sess1,'gradm',self.gradm)

    931. 

  931.         mlabraw.put(sess1,'hessm',self.hessm)

    932. 

  932. #----------------------------------------------------------------------------------------------------------------------

    933. 

  933. class MatWood:

    934. 

  934. 

    935. 

  935.     def __init__(self,intup):

    936. 

  936.         self.jAA = intup[0]

    937. 

  937.         self.jBB = intup[1]

    938. 

  938.         self.nexo = intup[2]

    939. 

  939.         self.ncon = intup[3]

    940. 

  940.         self.nendo = intup[4]

    941. 

  941.         self.sess1 = intup[5]

    942. 

  942.         self.NY = self.ncon+self.nendo

    943. 

  943.         self.NK = self.nendo

    944. 

  944.         self.NX = self.nexo

    945. 

  945.         self.mkwoodm()

    946. 

  946. 

    947. 

  947.     def solve(self):

    948. 

  948.         self.reds()

    949. 

  949.         self.solds()

    950. 

  950. 

    951. 

  951.     def mkwoodm(self):

    952. 

  952.         AA = self.jAA

    953. 

  953.         BB = self.jBB

    954. 

  954.         nendo = self.nendo

    955. 

  955.         nexo = self.nexo

    956. 

  956.         ncon = self.ncon

    957. 

  957.         nstates = nexo+nendo

    958. 

  958. 

    959. 

  959.         wAA = MAT.hstack((AA[:-nexo,nstates:],AA[:-nexo,nexo:nstates]))

    960. 

  960.         #wAA = MAT.hstack((AA[:,nstates:],AA[:,:nstates]))

    961. 

  961.         wBB = MAT.hstack((BB[:-nexo,nstates:],BB[:-nexo,nexo:nstates]))

    962. 

  962.         #wBB = MAT.hstack((BB[:,nstates:],BB[:,:nstates]))

    963. 

  963.         wCC = BB[:-nexo,:nexo]

    964. 

  964.         #wCC = MAT.zeros((ncon+nstates,1))

    965. 

  965. 

    966. 

  966. 

    967. 

  967.         self.wAA = wAA

    968. 

  968.         self.wBB = wBB

    969. 

  969.         self.wCC = wCC

    970. 

  970. 

    971. 

  971.     def reds(self):

    972. 

  972.         A = self.wAA

    973. 

  973.         B = self.wBB

    974. 

  974.         C = self.wCC

    975. 

  975.         NX = self.NX

    976. 

  976.         NK = self.NK

    977. 

  977.         NY = self.NY

    978. 

  978.         sess1 = self.sess1

    979. 

  979.         mlabraw.eval(sess1,'clear all;')

    980. 

  980.         mlabraw.eval(sess1,'cd '+mlabpath)

    981. 

  981.         mlabraw.eval(sess1,'cd Woodford')

    982. 

  982.         mlabraw.put(sess1,'wAA',A)

    983. 

  983.         mlabraw.put(sess1,'wBB',B)

    984. 

  984.         mlabraw.put(sess1,'wCC',C)

    985. 

  985.         mlabraw.put(sess1,'NX',NX)

    986. 

  986.         mlabraw.put(sess1,'NY',NY)

    987. 

  987.         mlabraw.put(sess1,'NK',NK)

    988. 

  988.         mlabraw.eval(sess1,'[Br,Cr,Lr,NF] = redsf(wAA,wBB,wCC,NY,NX,NK)')

    989. 

  989.         self.Br = np.matrix(mlabraw.get(sess1,'Br'))

    990. 

  990.         self.Cr = np.matrix(mlabraw.get(sess1,'Cr'))

    991. 

  991.         self.Lr = np.matrix(mlabraw.get(sess1,'Lr'))

    992. 

  992.         self.NF = np.matrix(mlabraw.get(sess1,'NF'))

    993. 

  993. 

    994. 

  994.     def solds(self):

    995. 

  995.         Br = self.Br

    996. 

  996.         Cr = self.Cr

    997. 

  997.         Lr = self.Lr

    998. 

  998.         NF = self.NF

    999. 

  999.         NX = self.NX

    1000. 

  1000.         NK = self.NK

    1001. 

  1001.         NY = self.NY

    1002. 

  1002.         sess1 = self.sess1

    1003. 

  1003.         mlabraw.eval(sess1,'clear all;')

    1004. 

  1004.         mlabraw.eval(sess1,'cd '+mlabpath)

    1005. 

  1005.         mlabraw.eval(sess1,'cd Woodford')

    1006. 

  1006.         mlabraw.put(sess1,'Br',Br)

    1007. 

  1007.         mlabraw.put(sess1,'Cr',Cr)

    1008. 

  1008.         mlabraw.put(sess1,'Lr',Lr)

    1009. 

  1009.         mlabraw.put(sess1,'NF',NF)

    1010. 

  1010.         mlabraw.put(sess1,'NX',NX)

    1011. 

  1011.         mlabraw.put(sess1,'NY',NY)

    1012. 

  1012.         mlabraw.put(sess1,'NK',NK)

    1013. 

  1013.         mlabraw.eval(sess1,'[D,F,G,H] = soldsf(Br,Cr,Lr,NY,NX,NK,NF)')

    1014. 

  1014.         self.D = np.matrix(mlabraw.get(sess1,'D'))

    1015. 

  1015.         self.F = np.matrix(mlabraw.get(sess1,'F'))

    1016. 

  1016.         self.G = np.matrix(mlabraw.get(sess1,'G'))

    1017. 

  1017.         self.H = np.matrix(mlabraw.get(sess1,'H'))

    1018. 

  1018. #----------------------------------------------------------------------------------------------------------------------

    1019. 

  1019. class ForKlein:

    1020. 

  1020. 

    1021. 

  1021.     def __init__(self,intup):

    1022. 

  1022.         self.uhlig = PyUhlig(intup)

    1023. 

  1023.         self.uhlig.mkmats()

    1024. 

  1024.         self.AA = self.uhlig.AA

    1025. 

  1025.         self.BB = self.uhlig.BB

    1026. 

  1026.         self.CC = self.uhlig.CC

    1027. 

  1027.         self.DD = self.uhlig.DD

    1028. 

  1028.         self.FF = self.uhlig.FF

    1029. 

  1029.         self.GG = self.uhlig.GG

    1030. 

  1030.         self.HH = self.uhlig.HH

    1031. 

  1031.         self.JJ = self.uhlig.JJ

    1032. 

  1032.         self.KK = self.uhlig.KK

    1033. 

  1033.         self.LL = self.uhlig.LL

    1034. 

  1034.         self.MM = self.uhlig.MM

    1035. 

  1035.         self.NN = self.uhlig.NN

    1036. 

  1036.         self.mkKleinMats()

    1037. 

  1037.         self.outdic = {}

    1038. 

  1038. 

    1039. 

  1039.     def mkKleinMats(self):

    1040. 

  1040.         # Make uhlig matrices locally available for computations

    1041. 

  1041.         AA = self.AA

    1042. 

  1042.         BB = self.BB

    1043. 

  1043.         CC = self.CC

    1044. 

  1044.         DD = self.DD

    1045. 

  1045.         FF = self.FF

    1046. 

  1046.         GG = self.GG

    1047. 

  1047.         HH = self.HH

    1048. 

  1048.         JJ = self.JJ

    1049. 

  1049.         KK = self.KK

    1050. 

  1050.         LL = self.LL

    1051. 

  1051.         MM = self.MM

    1052. 

  1052.         NN = self.NN

    1053. 

  1053. 

    1054. 

  1054.         # Determine size of states, endogenous

    1055. 

  1055.         exo_st = MAT.shape(NN)[1]

    1056. 

  1056.         endo_st = MAT.shape(BB)[1]

    1057. 

  1057.         endo_cn = MAT.shape(CC)[1]

    1058. 

  1058.         n_deteq = MAT.shape(AA)[0]

    1059. 

  1059.         n_expeq = MAT.shape(JJ)[0]

    1060. 

  1060.         tot_st = exo_st+endo_st

    1061. 

  1061.         self.tstates = tot_st

    1062. 

  1062.         tot_var = tot_st+endo_cn

    1063. 

  1063. 

    1064. 

  1064.         klein_A_rtwo = MAT.hstack((LL,GG))

    1065. 

  1065.         klein_A_rtwo = MAT.hstack((klein_A_rtwo,JJ))

    1066. 

  1066.         klein_A_rtwo_rows = MAT.shape(klein_A_rtwo)[0]

    1067. 

  1067.         klein_A_rtwo_cols = MAT.shape(klein_A_rtwo)[1]

    1068. 

  1068. 

    1069. 

  1069.         klein_A_rone = MAT.zeros((exo_st,klein_A_rtwo_cols))

    1070. 

  1070.         klein_A_rone = MAT.hstack((MAT.identity(exo_st),klein_A_rone[:,exo_st:]))

    1071. 

  1071. 

    1072. 

  1072.         klein_A_rthree = MAT.hstack((MAT.zeros((n_deteq,exo_st)),AA))

    1073. 

  1073.         klein_A_rthree = MAT.hstack((klein_A_rthree,MAT.zeros((n_deteq,endo_cn))))

    1074. 

  1074. 

    1075. 

  1075.         klein_A = MAT.vstack((klein_A_rone,klein_A_rtwo))

    1076. 

  1076.         klein_A = MAT.vstack((klein_A,klein_A_rthree))

    1077. 

  1077. 

    1078. 

  1078.         klein_B_rone = MAT.zeros((exo_st,klein_A_rtwo_cols))

    1079. 

  1079. 

    1080. 

  1080.         klein_B_rone = MAT.hstack((NN,klein_B_rone[:,exo_st:]))

    1081. 

  1081.         klein_B_rtwo = MAT.hstack((-MM,-HH))

    1082. 

  1082.         klein_B_rtwo = MAT.hstack((klein_B_rtwo,-KK))

    1083. 

  1083.         klein_B_rthree = MAT.hstack((-DD,-BB))

    1084. 

  1084.         klein_B_rthree = MAT.hstack((klein_B_rthree,-CC))

    1085. 

  1085. 

    1086. 

  1086.         klein_B = MAT.vstack((klein_B_rone,klein_B_rtwo))

    1087. 

  1087.         klein_B = MAT.vstack((klein_B,klein_B_rthree))  

    1088. 

  1088. 

    1089. 

  1089.         self.A = klein_A

    1090. 

  1090.         self.B = klein_B

    1091. 

  1091. 

    1092. 

  1092.     def solve(self):

    1093. 

  1093.         A = self.A

    1094. 

  1094.         B = self.B

    1095. 

  1095.         tstates = MAT.shape(self.AA)[1] + MAT.shape(self.DD)[1]

    1096. 

  1096.         (F,P,retcon) = isolab(A,B,tstates,MAT.shape(A)[0])

    1097. 

  1097.         if MAT.sum(P.reshape(-1,1)) == 0.0:

    1098. 

  1098.             return

    1099. 

  1099.         else:

    1100. 

  1100.             self.P = np.matrix(P)

    1101. 

  1101.             self.F = np.matrix(F)

    1102. 

  1102. #----------------------------------------------------------------------------------------------------------------------

    1103. 

  1103. class PyKlein2D(object):

    1104. 

  1104. 

    1105. 

  1105.     def __init__(self,intup,other=None):

    1106. 

  1106.         if other != None:

    1107. 

  1107.             self.other = other

    1108. 

  1108.         self.gra = intup[0]

    1109. 

  1109.         self.hes = intup[1]

    1110. 

  1110.         self.nendo = intup[2]

    1111. 

  1111.         self.nexo = intup[3]

    1112. 

  1112.         self.ncon = intup[4]

    1113. 

  1113.         self.sigma = intup[5]

    1114. 

  1114.         self.A = intup[6]

    1115. 

  1115.         self.B = intup[7]

    1116. 

  1116.         self.vardic = intup[8]

    1117. 

  1117.         self.vdic = intup[9]

    1118. 

  1118.         self.modname = intup[10]

    1119. 

  1119.         self.audic = intup[11]

    1120. 

  1120.         self.nacon = len(self.audic['con']['var'])

    1121. 

  1121.         self.naendo = len(self.audic['endo']['var'])

    1122. 

  1122.         if self.vardic['other']['var']:

    1123. 

  1123.             self.numjl = intup[12]

    1124. 

  1124.             self.numhl = intup[13]

    1125. 

  1125.             self.nother = intup[14]

    1126. 

  1126.             self.oswitch = 1

    1127. 

  1127.             kintup = (self.gra,self.nendo,

    1128. 

  1128.                   self.nexo,self.ncon,

    1129. 

  1129.                   self.sigma,self.A,self.B,

    1130. 

  1130.                   self.vardic,self.vdic,

    1131. 

  1131.                   self.modname,self.audic,

    1132. 

  1132.                   self.numjl,self.nother)

    1133. 

  1133.         else:

    1134. 

  1134.             self.oswitch = 0

    1135. 

  1135.             kintup = (self.gra,self.nendo,

    1136. 

  1136.                   self.nexo,self.ncon,

    1137. 

  1137.                   self.sigma,self.A,self.B,

    1138. 

  1138.                   self.vardic,self.vdic,

    1139. 

  1139.                   self.modname,self.audic)

    1140. 

  1140.         self.tstates = self.nendo+self.nexo

    1141. 

  1141.         self.forkleind = ForKleinD(kintup)

    1142. 

  1142.         self.ssigma = self.mkssigma(self.tstates,self.nexo,self.sigma)

    1143. 

  1143. 

    1144. 

  1144.     def mkssigma(self,tstates,nexo,sigma):

    1145. 

  1145.         ssigma = MAT.zeros((tstates,tstates))

    1146. 

  1146.         for i in range(nexo):

    1147. 

  1147.             ssigma[i,i] = sigma[i,i]

    1148. 

  1148.         return ssigma

    1149. 

  1149. 

    1150. 

  1150.     def solab2(self):

    1151. 

  1151. 

    1152. 

  1152.         # Tracem function

    1153. 

  1153.         def tracem(xmat):

    1154. 

  1154.             n = MAT.shape(xmat)[1]

    1155. 

  1155.             m = MAT.shape(xmat)[0]/n

    1156. 

  1156.             ymat = MAT.zeros((m,1))

    1157. 

  1157.             for i1 in xrange(int(m)):

    1158. 

  1158.                 ymat[i1,0]=MAT.trace(xmat[n*i1:i1*n+1,:n])

    1159. 

  1159.             return ymat

    1160. 

  1160. 

    1161. 

  1161.         ssigma = self.ssigma

    1162. 

  1162.         pp = self.P

    1163. 

  1163.         ff = self.F

    1164. 

  1164.         gra = self.gra

    1165. 

  1165.         hes = self.hes

    1166. 

  1166.         m = self.ncon+self.nendo+self.nexo

    1167. 

  1167.         ny = self.ncon

    1168. 

  1168.         nx = self.tstates

    1169. 

  1169. 

    1170. 

  1170.         f1 = gra[:,:nx]

    1171. 

  1171.         f2 = gra[:,nx:m]

    1172. 

  1172.         f4 = gra[:,m+nx:2*m]

    1173. 

  1173. 

    1174. 

  1174.         mm = MAT.vstack((pp,ff*pp,MAT.eye(nx),ff))

    1175. 

  1175.         aa1 = MAT.kron(MAT.eye(m),mm.T)*hes*mm

    1176. 

  1176.         bb1 = MAT.kron(f1,MAT.eye(nx))

    1177. 

  1177.         bb2 = MAT.kron(f2,MAT.eye(nx))

    1178. 

  1178.         bb4 = MAT.kron(f4,MAT.eye(nx))

    1179. 

  1179.         cc1 = MAT.kron(MAT.eye(ny),pp.T)

    1180. 

  1180.         cc2 = MAT.kron(ff,MAT.eye(nx))

    1181. 

  1181. 

    1182. 

  1182.         aa_1 = MAT.kron(MAT.eye(nx),bb4)+MAT.kron(pp.T,bb2*cc1)

    1183. 

  1183.         aa_2 = MAT.kron(MAT.eye(nx),bb1+bb2*cc2)

    1184. 

  1184.         aa = MAT.hstack((aa_1,aa_2))

    1185. 

  1185.         sol = np.linalg.solve(-aa,HLP.ctr(aa1.flatten(1)))

    1186. 

  1186.         ee = HLP.ctr(sol[:nx**2*ny].reshape(-1,nx*ny))

    1187. 

  1187.         gg = HLP.ctr(sol[nx**2*ny:].reshape(-1,nx**2))

    1188. 

  1188.         ma = 2.0*MAT.hstack((f1+f2*ff,f2+f4))

    1189. 

  1189.         eyeff = MAT.vstack((MAT.eye(nx),ff,MAT.zeros((m,nx))))

    1190. 

  1190.         ve_1 = f2*tracem(MAT.kron(MAT.eye(ny),ssigma)*ee)

    1191. 

  1191.         ve_2 = tracem(MAT.kron(MAT.eye(m),HLP.ctr(eyeff))*hes*eyeff*ssigma)

    1192. 

  1192.         ve = ve_1 + ve_2

    1193. 

  1193.         kxy = -(ma.I)*ve

    1194. 

  1194.         kx = kxy[:nx]

    1195. 

  1195.         ky = kxy[nx:]

    1196. 

  1196.         self.E = ee

    1197. 

  1197.         self.G = gg

    1198. 

  1198.         self.KX = kx

    1199. 

  1199.         self.KY = ky

    1200. 

  1200. 

    1201. 

  1201.     def solve(self):

    1202. 

  1202.         self.forkleind.solve()

    1203. 

  1203.         self.P = self.forkleind.P

    1204. 

  1204.         self.F = self.forkleind.F

    1205. 

  1205.         self.solab2()

    1206. 

  1206. 

    1207. 

  1207.     def sim(self,tlen,sntup=None):

    1208. 

  1208.         # Deals with 0-tuples

    1209. 

  1209.         if type(sntup) != type((1,2,3)) and type(sntup) == type('abc'):

    1210. 

  1210.             sntup = (sntup,)

    1211. 

  1211.         # Add 1000 observations, as they will be thrown away

    1212. 

  1212.         # Add one more observation to start first-order vector

    1213. 

  1213.         exoli = [x[1] for x in self.vardic['exo']['var']]

    1214. 

  1214.         indx=[]

    1215. 

  1215.         if sntup == None:

    1216. 

  1216.             indx = range(len(exoli))

    1217. 

  1217.         else:

    1218. 

  1218.             for name in sntup:

    1219. 

  1219.                 if name not in exoli:

    1220. 

  1220.                     return 'Error: '+name+' is not a valid exo variable for this model!'

    1221. 

  1221.                 else:

    1222. 

  1222.                     indx.append(exoli.index(name))

    1223. 

  1223.         indx.sort()

    1224. 

  1224.         ncon = self.ncon

    1225. 

  1225.         nexo = self.nexo

    1226. 

  1226.         nendo = self.nendo

    1227. 

  1227.         tstates = self.tstates

    1228. 

  1228.         tlena = 1000+tlen

    1229. 

  1229.         sigma = self.sigma

    1230. 

  1230.         ssigma = self.ssigma

    1231. 

  1231.         if self.oswitch:

    1232. 

  1232.             numjl = self.numjl

    1233. 

  1233.             numhl = self.numhl

    1234. 

  1234.         kx = self.KX

    1235. 

  1235.         ky = self.KY

    1236. 

  1236.         pp = self.P

    1237. 

  1237.         ff = self.F

    1238. 

  1238.         ee = self.E

    1239. 

  1239.         gg = self.G

    1240. 

  1240.         count = 0

    1241. 

  1241.         for varia in MAT.diag(sigma):

    1242. 

  1242.             if locals().has_key('ranvec'):

    1243. 

  1243.                 if count in indx:

    1244. 

  1244.                     ranvec = MAT.vstack((ranvec,N.sqrt(varia)*MAT.matrix(N.random.standard_normal(tlena))))

    1245. 

  1245.                 else:

    1246. 

  1246.                     ranvec = MAT.vstack((ranvec,MAT.zeros((1,tlena))))

    1247. 

  1247.             else:

    1248. 

  1248.                 if count in indx:

    1249. 

  1249.                     ranvec = np.sqrt(varia)*MAT.matrix(N.random.standard_normal(tlena))

    1250. 

  1250.                 else:

    1251. 

  1251.                     ranvec = MAT.zeros((1,tlena))

    1252. 

  1252.             count = count + 1

    1253. 

  1253. 

    1254. 

  1254.         ranvec = MAT.vstack((ranvec,MAT.zeros((nendo,tlena))))

    1255. 

  1255. 

    1256. 

  1256.         x_one_m1 = ranvec[:,0]

    1257. 

  1257.         x_one_0 = pp*x_one_m1

    1258. 

  1258.         x_two_m1 = kx + ranvec[:,0]

    1259. 

  1259.         y_one_0 = ff*x_one_m1

    1260. 

  1260.         y_one_m1 = MAT.zeros(y_one_0.shape)

    1261. 

  1261.         y_two_0 = ky+0.5*MAT.kron(MAT.eye(ncon),x_one_m1.T)*ee*x_one_m1

    1262. 

  1262.         if self.oswitch:

    1263. 

  1263.             nother = self.nother

    1264. 

  1264.             o_one_0 = numjl*MAT.vstack((x_one_0,y_one_0,x_one_m1,y_one_m1))

    1265. 

  1265.             o_two_0 = numjl*MAT.vstack((x_one_0,y_one_0,x_one_m1,y_one_m1))+\

    1266. 

  1266.                 0.5*MAT.kron(MAT.eye(nother),MAT.vstack((x_one_0,y_one_0,x_one_m1,y_one_m1)).T)*\

    1267. 

  1267.                 numhl*MAT.vstack((x_one_0,y_one_0,x_one_m1,y_one_m1))

    1268. 

  1268. 

    1269. 

  1269.         x_one_c = COP.deepcopy(x_one_m1)

    1270. 

  1270.         y_one_c = COP.deepcopy(y_one_0)

    1271. 

  1271.         x_two_c = COP.deepcopy(x_two_m1)

    1272. 

  1272.         y_two_c = COP.deepcopy(y_two_0)

    1273. 

  1273.         if self.oswitch:

    1274. 

  1274.             o_one_c = COP.deepcopy(o_one_0)

    1275. 

  1275.             o_two_c = COP.deepcopy(o_two_0)

    1276. 

  1276. 

    1277. 

  1277.         x_one = COP.deepcopy(x_one_c)

    1278. 

  1278.         y_one = COP.deepcopy(y_one_c)

    1279. 

  1279.         x_two = COP.deepcopy(x_two_c)

    1280. 

  1280.         y_two = COP.deepcopy(y_two_c)

    1281. 

  1281.         if self.oswitch:

    1282. 

  1282.             o_one = COP.deepcopy(o_one_c)

    1283. 

  1283.             o_two = COP.deepcopy(o_two_c)

    1284. 

  1284. 

    1285. 

  1285.         for i1 in range(1,ranvec.shape[1],1):

    1286. 

  1286. 

    1287. 

  1287.             x_one_n = pp*x_one_c+ranvec[:,i1]

    1288. 

  1288.             y_one_n = ff*x_one_c

    1289. 

  1289. 

    1290. 

  1290.             x_two_n = kx+pp*x_two_c+\

    1291. 

  1291.                 0.5*MAT.kron(MAT.eye(tstates),x_one_c.T)*\

    1292. 

  1292.                 gg*x_one_c+ranvec[:,i1]

    1293. 

  1293. 

    1294. 

  1294.             y_two_n = ky+ff*x_two_c+\

    1295. 

  1295.                 0.5*MAT.kron(MAT.eye(ncon),x_one_c.T)*\

    1296. 

  1296.                 ee*x_one_c

    1297. 

  1297. 

    1298. 

  1298.             if self.oswitch:

    1299. 

  1299.                 o_one_n = numjl*MAT.vstack((x_one_n,y_one_n,x_one_c,y_one_c))

    1300. 

  1300.                 o_two_n = numjl*MAT.vstack((x_one_n,y_one_n,x_one_c,y_one_c))+\

    1301. 

  1301.                     0.5*MAT.kron(MAT.eye(nother),MAT.vstack((x_one_n,y_one_n,x_one_c,y_one_c)).T)*\

    1302. 

  1302.                     numhl*MAT.vstack((x_one_n,y_one_n,x_one_c,y_one_c))

    1303. 

  1303. 

    1304. 

  1304.             x_one = MAT.hstack((x_one,x_one_c))

    1305. 

  1305.             y_one = MAT.hstack((y_one,y_one_n))

    1306. 

  1306.             x_two = MAT.hstack((x_two,x_two_c))

    1307. 

  1307.             y_two = MAT.hstack((y_two,y_two_n))

    1308. 

  1308.             if self.oswitch:

    1309. 

  1309.                 o_one = MAT.hstack((o_one,o_one_n))

    1310. 

  1310.                 o_two = MAT.hstack((o_two,o_two_n))

    1311. 

  1311. 

    1312. 

  1312.             x_one_c = x_one_n

    1313. 

  1313.             y_one_c = y_one_n

    1314. 

  1314.             x_two_c = x_one_n

    1315. 

  1315.             y_two_c = y_two_n

    1316. 

  1316.             if self.oswitch:

    1317. 

  1317.                 o_one_c = o_one_n

    1318. 

  1318.                 o_two_c = o_two_n

    1319. 

  1319. 

    1320. 

  1320.         # Throw away first 1000 obs

    1321. 

  1321.         x_one = x_one[:,1000:]

    1322. 

  1322.         y_one = y_one[:,1000:]

    1323. 

  1323.         y_two = y_two[:,1000:]

    1324. 

  1324.         x_two = x_two[:,1000:]

    1325. 

  1325.         if self.oswitch:

    1326. 

  1326.             o_one = o_one[:,1000:]

    1327. 

  1327.             o_two = o_two[:,1000:]

    1328. 

  1328. 

    1329. 

  1329.         self.sim_y_two = y_two

    1330. 

  1330.         self.sim_y_one = y_one

    1331. 

  1331.         self.sim_x_two = x_two

    1332. 

  1332.         self.sim_x_one = x_one

    1333. 

  1333.         self.insim = [self.sim_x_two,self.sim_y_two]

    1334. 

  1334.         if self.oswitch:

    1335. 

  1335.             self.sim_o_one = o_one

    1336. 

  1336.             self.sim_o_two = o_two

    1337. 

  1337.             self.insim = self.insim + [self.sim_o_two,]

    1338. 

  1338. 

    1339. 

  1339.     def mkcvar(self,vname=False,insim='insim'):

    1340. 

  1340.         # Now also produce unconditional variances, and normal. variances with v

    1341. 

  1341.         # Check if simul data is available

    1342. 

  1342.         if insim not in dir(self):

    1343. 

  1343.             return "Error: Can't produce uncond. variance table without simulation data!"

    1344. 

  1344. 

    1345. 

  1345.         exoli = [x[1] for x in self.vardic['exo']['var']]

    1346. 

  1346.         endoli = [x[1] for x in self.vardic['endo']['var']]

    1347. 

  1347.         stateli = exoli + endoli

    1348. 

  1348.         conli = [x[1] for x in self.vardic['con']['var']]

    1349. 

  1349.         alli = stateli+conli

    1350. 

  1350.         if self.oswitch:

    1351. 

  1351.             otherli = [x[1] for x in self.vardic['other']['var']]

    1352. 

  1352.             alli = alli + otherli

    1353. 

  1353. 

    1354. 

  1354.         # Check if vname is valid

    1355. 

  1355.         if vname not in alli:

    1356. 

  1356.             return 'Error: '+vname+' is not a valid variable for this model!'

    1357. 

  1357.         insim = eval('self.'+insim)

    1358. 

  1358.         osim_x = COP.deepcopy(insim[0])

    1359. 

  1359.         osim_y = COP.deepcopy(insim[1])

    1360. 

  1360.         nendo = self.nendo

    1361. 

  1361.         nexo = self.nexo

    1362. 

  1362.         ncon = self.ncon

    1363. 

  1363. 

    1364. 

  1364.         # Now hp filter the simulations before graphing according to filtup

    1365. 

  1365.         for i1 in xrange(osim_y.shape[0]):

    1366. 

  1366.             yy = osim_y[i1,:].__array__().T

    1367. 

  1367.   …
    1368.

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