PageRenderTime 63ms CodeModel.GetById 17ms RepoModel.GetById 0ms app.codeStats 1ms

/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 
    

  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.             woy = np.zeros((osim_y.shape[1],3))

    1368. 

  1368.             lam = 1600

    1369. 

  1369.             yyf = MAT.matrix(hpfilter(data=yy,lam=1600))

    1370. 

  1370.             osim_y[i1,:] = yyf[0]

    1371. 

  1371.         # Now filter the state variables!

    1372. 

  1372.         for i1 in xrange(osim_x.shape[0]):

    1373. 

  1373.             xx = osim_x[i1,:].__array__().T

    1374. 

  1374.             wox = np.zeros((osim_x.shape[1],3))

    1375. 

  1375.             lam = 1600

    1376. 

  1376.             xxf = MAT.matrix(hpfilter(data=xx,lam=1600))

    1377. 

  1377.             osim_x[i1,:] = xxf[0]

    1378. 

  1378. 

    1379. 

  1379.         if self.oswitch:

    1380. 

  1380.             osim_o = COP.deepcopy(insim[2])

    1381. 

  1381.             # Now hp filter the other variables!

    1382. 

  1382.             for i1 in xrange(osim_o.shape[0]):

    1383. 

  1383.                 oo = osim_o[i1,:].__array__().T

    1384. 

  1384.                 woo = np.zeros((osim_o.shape[1],3))

    1385. 

  1385.                 lam = 1600

    1386. 

  1386.                 oof = MAT.matrix(hpfilter(data=oo,lam=1600))

    1387. 

  1387.                 osim_o[i1,:] = oof[0]

    1388. 

  1388. 

    1389. 

  1389.         # Stack all matrices into one

    1390. 

  1390.         allmat = MAT.vstack((osim_x,osim_y))

    1391. 

  1391.         if self.oswitch:

    1392. 

  1392.             allmat = MAT.vstack((allmat,osim_o))

    1393. 

  1393.         varmat = MAT.zeros((allmat.shape[0],1))

    1394. 

  1394.         for i1 in range(allmat.shape[0]):

    1395. 

  1395.             varmat[i1,0] = np.var(allmat[i1,:].A.flatten(1))

    1396. 

  1396. 

    1397. 

  1397.         if vname:

    1398. 

  1398.             relmat = MAT.zeros((allmat.shape[0],1))

    1399. 

  1399.             indx = alli.index(vname)

    1400. 

  1400.             varvari = np.var(allmat[indx,:].A.flatten(1))

    1401. 

  1401.             for i1 in range(allmat.shape[0]):

    1402. 

  1402.                 relmat[i1,0] = varmat[i1,0]/varvari

    1403. 

  1403.             self.cvarm = MAT.hstack((varmat,relmat))

    1404. 

  1404.         else:

    1405. 

  1405.             self.cvarm = varmat

    1406. 

  1406. 

    1407. 

  1407. 

    1408. 

  1408. 

    1409. 

  1409.     def mkact(self,vname,intup=None,insim='insim'):

    1410. 

  1410.         # Now also produce autocorrelation table

    1411. 

  1411.         # Check if simul data is available

    1412. 

  1412.         if insim not in dir(self):

    1413. 

  1413.             return "Error: Can't produce autocorrelation table without simulation data!"

    1414. 

  1414. 

    1415. 

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

    1416. 

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

    1417. 

  1417.         stateli = exoli + endoli

    1418. 

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

    1419. 

  1419.         alli = stateli+conli

    1420. 

  1420.         if self.oswitch:

    1421. 

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

    1422. 

  1422.             alli = alli + otherli

    1423. 

  1423.         # Check if vname is valid

    1424. 

  1424.         if vname not in alli:

    1425. 

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

    1426. 

  1426. 

    1427. 

  1427.         if intup == None:

    1428. 

  1428.             lags = 3

    1429. 

  1429.             leads = 3

    1430. 

  1430.         else:

    1431. 

  1431.             lags = intup[0]

    1432. 

  1432.             leads = intup[1]

    1433. 

  1433. 

    1434. 

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

    1435. 

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

    1436. 

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

    1437. 

  1437. 

    1438. 

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

    1439. 

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

    1440. 

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

    1441. 

  1441.             woy = np.zeros((osim_y.shape[1],3))

    1442. 

  1442.             lam = 1600

    1443. 

  1443.             yyf = MAT.matrix(hpfilter(data=yy,lam=1600))

    1444. 

  1444.             osim_y[i1,:] = yyf[0]

    1445. 

  1445.         # Now filter the state variables!

    1446. 

  1446.         for i1 in xrange(osim_x.shape[0]):

    1447. 

  1447.             xx = osim_x[i1,:].__array__().T

    1448. 

  1448.             wox = np.zeros((osim_x.shape[1],3))

    1449. 

  1449.             lam = 1600

    1450. 

  1450.             xxf = MAT.matrix(hpfilter(data=xx,lam=1600))

    1451. 

  1451.             osim_x[i1,:] = xxf[0]

    1452. 

  1452. 

    1453. 

  1453.         sim_x = osim_x[:,lags:-leads]

    1454. 

  1454.         sim_xf = osim_x[:,leads+1:]

    1455. 

  1455.         sim_yf = osim_y[:,leads+1:]

    1456. 

  1456.         sim_y = osim_y[:,lags:-leads]

    1457. 

  1457.         sim_xl = osim_x[:,:-lags-1]

    1458. 

  1458.         sim_yl = osim_y[:,:-lags-1]

    1459. 

  1459. 

    1460. 

  1460.         if self.oswitch:

    1461. 

  1461.             osim_o = COP.deepcopy(insim[2])

    1462. 

  1462.             # Now hp filter the other variables!

    1463. 

  1463.             for i1 in xrange(osim_o.shape[0]):

    1464. 

  1464.                 oo = osim_o[i1,:].__array__().T

    1465. 

  1465.                 woo = np.zeros((osim_o.shape[1],3))

    1466. 

  1466.                 lam = 1600

    1467. 

  1467.                 oof = MAT.matrix(hpfilter(data=oo,lam=1600))

    1468. 

  1468.                 osim_o[i1,:] = oof[0]

    1469. 

  1469.             sim_o = osim_o[:,lags:-leads]

    1470. 

  1470.             sim_of = osim_o[:,leads+1:]

    1471. 

  1471.             sim_ol = osim_o[:,:-lags-1]

    1472. 

  1472.             alli = alli + otherli

    1473. 

  1473. 

    1474. 

  1474.         posa = alli.index(vname)

    1475. 

  1475.         actm = MAT.zeros((len(alli),lags+1+leads))

    1476. 

  1476.         sima = MAT.vstack((sim_x,sim_y))

    1477. 

  1477.         sima_f = MAT.vstack((sim_xf,sim_yf))

    1478. 

  1478.         sima_l = MAT.vstack((sim_xl,sim_yl))

    1479. 

  1479.         if self.oswitch:

    1480. 

  1480.             sima = MAT.vstack((sima,sim_o))

    1481. 

  1481.             sima_f = MAT.vstack((sima_f,sim_of))

    1482. 

  1482.             sima_l = MAT.vstack((sima_l,sim_ol))

    1483. 

  1483.         for i1 in range(sima.shape[0]):

    1484. 

  1484.             # Do lags

    1485. 

  1485.             for i2 in range(lags):

    1486. 

  1486.                 actm[i1,i2] = np.round(np.corrcoef(sima[posa,:].A.flatten(1),sima_l[i1,i2:sima_l.shape[1]-lags+i2+1].A.flatten(1))[1,0],3)

    1487. 

  1487.             # Do current

    1488. 

  1488.             actm[i1,lags] = np.round(np.corrcoef(sima[posa,:].A.flatten(1),sima[i1,:].A.flatten(1))[1,0],3)

    1489. 

  1489.             # Do leads

    1490. 

  1490.             for i2 in range(leads):

    1491. 

  1491.                 actm[i1,lags+1+i2] = np.round(np.corrcoef(sima[posa,:].A.flatten(1),sima_f[i1,i2:sima_f.shape[1]-leads+i2+1].A.flatten(1))[1,0],3)

    1492. 

  1492. 

    1493. 

  1493.         self.actm = actm

    1494. 

  1494.         self.actname = vname

    1495. 

  1495.         if intup:

    1496. 

  1496.             self.actin = intup

    1497. 

  1497. 

    1498. 

  1498.     def show_act(self):

    1499. 

  1499.         if 'actname' not in dir(self):

    1500. 

  1500.             return 'Error: You have not produced any autocorrelation table yet!'

    1501. 

  1501.         nexo = self.nexo

    1502. 

  1502.         nendo = self.nendo

    1503. 

  1503.         ncon = self.ncon

    1504. 

  1504.         if 'naendo' in dir(self):

    1505. 

  1505.             naendo = self.naendo

    1506. 

  1506.         else:

    1507. 

  1507.             naendo = 0

    1508. 

  1508.         if 'nacon' in dir(self):

    1509. 

  1509.             nacon = self.nacon

    1510. 

  1510.         else:

    1511. 

  1511.             nacon = 0

    1512. 

  1512.         nother = len(self.vardic['other']['var'])

    1513. 

  1513.         respva = COP.deepcopy(self.actname)

    1514. 

  1514.         if 'actin' in dir(self):

    1515. 

  1515.             lags = COP.deepcopy(int(self.actin[0]))

    1516. 

  1516.             leads = COP.deepcopy(int(self.actin[1]))

    1517. 

  1517.         else:

    1518. 

  1518.             lags = 3

    1519. 

  1519.             leads = 3

    1520. 

  1520.         actm = COP.deepcopy(self.actm)

    1521. 

  1521.         vdic = self.vdic

    1522. 

  1522.         allvari = [x[1] for x in vdic['exo']]

    1523. 

  1523.         allvari = allvari + [x[1] for x in vdic['endo']]

    1524. 

  1524.         allvari = allvari + [x[1] for x in vdic['con']]

    1525. 

  1525.         if self.oswitch:

    1526. 

  1526.             allvari = allvari + [x[1] for x in vdic['other']]

    1527. 

  1527.         allvari = allvari[:nexo+nendo-naendo]+allvari[nexo+nendo:nexo+nendo+ncon-nacon]

    1528. 

  1528.         actm = MAT.vstack((actm[:nexo+nendo-naendo,:],actm[nexo+nendo:nexo+nendo+ncon-nacon,:]))

    1529. 

  1529. 

    1530. 

  1530.         # determine longest varname

    1531. 

  1531.         vnmax = 0

    1532. 

  1532.         for vari in allvari:

    1533. 

  1533.             if len(vari) > vnmax:

    1534. 

  1534.                 vnmax = len(vari)

    1535. 

  1535.         print ''

    1536. 

  1536.         print 'Autocorrelation table, current '+respva

    1537. 

  1537.         print '='*65

    1538. 

  1538.         for i1,vari in enumerate(allvari):

    1539. 

  1539.             print vari+(vnmax-len(vari))*' '+'  |'+str(actm[i1,:])[2:-2]

    1540. 

  1540. 

    1541. 

  1541.     def show_sim(self,intup,filtup=None,insim='insim'):

    1542. 

  1542.         # Deals with 0-tuples

    1543. 

  1543.         if type(intup) != type((1,2,3)) and type(intup) == type('abc'):

    1544. 

  1544.             intup = (intup,)

    1545. 

  1545.         # Check if simulations have been carried out

    1546. 

  1546.         if insim not in dir(self):

    1547. 

  1547.             return 'Error: You have not produced any simulations yet! Nothing to graph!'

    1548. 

  1548.         mname = self.modname

    1549. 

  1549.         vardic = self.vardic

    1550. 

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

    1551. 

  1551.         sim_x = COP.deepcopy(insim[0])

    1552. 

  1552.         sim_y = COP.deepcopy(insim[1])

    1553. 

  1553.         tlen = sim_x.shape[1]

    1554. 

  1554.         if self.oswitch:

    1555. 

  1555.             sim_o = COP.deepcopy(insim[2])

    1556. 

  1556.             nother = self.nother

    1557. 

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

    1558. 

  1558.         nexo = self.nexo

    1559. 

  1559.         nendo = self.nendo

    1560. 

  1560.         ncon = self.ncon

    1561. 

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

    1562. 

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

    1563. 

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

    1564. 

  1564. 

    1565. 

  1565.         # If filtup is None, create filtup from vardic

    1566. 

  1566.         if filtup == None:

    1567. 

  1567.             filli = [0]*len(intup)

    1568. 

  1568.             for i1,x in enumerate(intup):

    1569. 

  1569.                 for y in vardic.keys():

    1570. 

  1570.                     if x in [z[1] for z in vardic[y]['var']]:

    1571. 

  1571.                         indx = [z[1] for z in vardic[y]['var']].index(x)

    1572. 

  1572.                         if 'hp' in [z for z in vardic[y]['mod']][indx]:

    1573. 

  1573.                             filli[i1] = 1

    1574. 

  1574.                         elif 'bk' in [z for z in vardic[y]['mod']][indx]:

    1575. 

  1575.                             filli[i1] = 2

    1576. 

  1576.                         elif 'cf' in [z for z in vardic[y]['mod']][indx]:

    1577. 

  1577.                             filli[i1] = 3

    1578. 

  1578.             filtup = tuple(filli)

    1579. 

  1579. 

    1580. 

  1580. 

    1581. 

  1581.         stateli = exoli+endoli

    1582. 

  1582.         alli = stateli+conli

    1583. 

  1583.         if self.oswitch:

    1584. 

  1584.             alli = alli+otherli

    1585. 

  1585.         # Check if all name in intup are available to graph

    1586. 

  1586.         for name in intup:

    1587. 

  1587.             if name not in alli:

    1588. 

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

    1589. 

  1589.         # Create x, y and z indeces

    1590. 

  1590.         indx = []

    1591. 

  1591.         indy = []

    1592. 

  1592.         indo = []

    1593. 

  1593.         for name in intup:

    1594. 

  1594.             if name in stateli:

    1595. 

  1595.                 indx.append(stateli.index(name))

    1596. 

  1596.             elif name in conli:

    1597. 

  1597.                 indy.append(conli.index(name))

    1598. 

  1598.             elif self.oswitch and name in otherli:

    1599. 

  1599.                 indo.append(otherli.index(name))

    1600. 

  1600.         leg = []

    1601. 

  1601.         indx.sort()

    1602. 

  1602.         indy.sort()

    1603. 

  1603.         indo.sort()     

    1604. 

  1604. 

    1605. 

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

    1606. 

  1606.         for i1 in xrange(sim_y.shape[0]):

    1607. 

  1607.             if indy and (i1 in indy) and filtup[list(intup).index(conli[i1])]:

    1608. 

  1608.                 if filtup[list(intup).index(conli[i1])] == 1:

    1609. 

  1609.                     conli[i1] = conli[i1]+'(hpf)'

    1610. 

  1610.                     yy = sim_y[i1,:].__array__().T

    1611. 

  1611.                     woy = np.zeros((tlen,3))

    1612. 

  1612.                     lam = 1600

    1613. 

  1613.                     yyf = MAT.matrix(hpfilter(data=yy,lam=1600))

    1614. 

  1614.                     sim_y[i1,:] = yyf[0]*100.0

    1615. 

  1615.                 elif filtup[list(intup).index(conli[i1])] == 2:

    1616. 

  1616.                     conli[i1] = conli[i1]+'(bkf)'

    1617. 

  1617.                     yy = sim_y[i1,:].__array__().T

    1618. 

  1618.                     woy = np.zeros((tlen,1))

    1619. 

  1619.                     up = 6

    1620. 

  1620.                     dn = 32

    1621. 

  1621.                     kkl = 12

    1622. 

  1622.                     yyf = MAT.matrix(bkfilter(data=yy,up=up,dn=dn,kkl=kkl))

    1623. 

  1623.                     sim_y[i1,:] = yyf[0]*100.0

    1624. 

  1624.                 elif filtup[list(intup).index(conli[i1])] == 3:

    1625. 

  1625.                     conli[i1] = conli[i1]+'(cff)'

    1626. 

  1626.                     yy = sim_y[i1,:].__array__().T

    1627. 

  1627.                     woy = np.zeros((tlen,1))

    1628. 

  1628.                     low = 6

    1629. 

  1629.                     high = 32

    1630. 

  1630.                     drift = True

    1631. 

  1631.                     yyf = MAT.matrix(cffilter(data=yy,low=low,high=high,drift=drift))

    1632. 

  1632.                     sim_y[i1,:] = yyf[0]*100.0

    1633. 

  1633.         # Now filter the state variables!

    1634. 

  1634.         for i1 in xrange(sim_x.shape[0]):

    1635. 

  1635.             if indx and (i1 in indx) and filtup[list(intup).index(stateli[i1])]:

    1636. 

  1636.                 if filtup[list(intup).index(stateli[i1])] == 1:

    1637. 

  1637.                     stateli[i1] = stateli[i1]+'(hpf)'

    1638. 

  1638.                     xx = sim_x[i1,:].__array__().T

    1639. 

  1639.                     wox = np.zeros((tlen,3))

    1640. 

  1640.                     lam = 1600

    1641. 

  1641.                     xxf = MAT.matrix(hpfilter(data=xx,lam=1600))

    1642. 

  1642.                     sim_x[i1,:] = xxf[0]*100.0

    1643. 

  1643.                 elif filtup[list(intup).index(stateli[i1])] == 2:

    1644. 

  1644.                     stateli[i1] = stateli[i1]+'(bkf)'

    1645. 

  1645.                     xx = sim_x[i1,:].__array__().T

    1646. 

  1646.                     wox = np.zeros((tlen,1))

    1647. 

  1647.                     up = 6

    1648. 

  1648.                     dn = 32

    1649. 

  1649.                     kkl = 12

    1650. 

  1650.                     xxf = MAT.matrix(bkfilter(data=xx,up=up,dn=dn,kkl=kkl))

    1651. 

  1651.                     sim_x[i1,:] = xxf[0]*100.0

    1652. 

  1652.                 elif filtup[list(intup).index(stateli[i1])] == 3:

    1653. 

  1653.                     stateli[i1] = stateli[i1]+'(cff)'

    1654. 

  1654.                     xx = sim_x[i1,:].__array__().T

    1655. 

  1655.                     wox = np.zeros((tlen,1))

    1656. 

  1656.                     low = 6

    1657. 

  1657.                     high = 32

    1658. 

  1658.                     drift = True

    1659. 

  1659.                     xxf = MAT.matrix(cffilter(data=xx,low=low,high=high,drift=drift))

    1660. 

  1660.                     sim_x[i1,:] = xxf[0]*100.0 

    1661. 

  1661.         # Now hp filter the other variables!

    1662. 

  1662.         if self.oswitch:

    1663. 

  1663.             for i1 in xrange(sim_o.shape[0]):

    1664. 

  1664.                 if indo and (i1 in indo) and filtup[list(intup).index(otherli[i1])]:

    1665. 

  1665.                     if filtup[list(intup).index(otherli[i1])] == 1:                   

    1666. 

  1666.                         otherli[i1] = otherli[i1]+'(hpf)'

    1667. 

  1667.                         oo = sim_o[i1,:].__array__().T

    1668. 

  1668.                         woo = np.zeros((tlen,3))

    1669. 

  1669.                         lam = 1600

    1670. 

  1670.                         oof = MAT.matrix(hpfilter(data=oo,lam=1600))

    1671. 

  1671.                     elif filtup[list(intup).index(otherli[i1])] == 2:

    1672. 

  1672.                         otherli[i1] = otherli[i1]+'(bkf)'

    1673. 

  1673.                         oo = sim_o[i1,:].__array__().T

    1674. 

  1674.                         woo = np.zeros((tlen,1))

    1675. 

  1675.                         up = 6

    1676. 

  1676.                         dn = 32

    1677. 

  1677.                         kkl = 12

    1678. 

  1678.                         oof = MAT.matrix(bkfilter(data=oo,up=up,dn=dn,kkl=kkl))

    1679. 

  1679.                     elif filtup[list(intup).index(otherli[i1])] == 3:

    1680. 

  1680.                         otherli[i1] = otherli[i1]+'(cff)'

    1681. 

  1681.                         oo = sim_o[i1,:].__array__().T

    1682. 

  1682.                         woo = np.zeros((tlen,1))

    1683. 

  1683.                         low = 6

    1684. 

  1684.                         high = 32

    1685. 

  1685.                         drift = True

    1686. 

  1686.                         oof = MAT.matrix(cffilter(data=oo,low=low,high=high,drift=drift)) 

    1687. 

  1687. 

    1688. 

  1688.         if indx and indy and indo:

    1689. 

  1689.             for x in indx:

    1690. 

  1690.                 leg.append(stateli[x])

    1691. 

  1691.             for y in indy:

    1692. 

  1692.                 leg.append(conli[y])

    1693. 

  1693.             for o in indo:

    1694. 

  1694.                 leg.append(otherli[o])

    1695. 

  1695.             leg = tuple(leg)

    1696. 

  1696.             figo = P.figure()

    1697. 

  1697.             P.plot(MAT.hstack((sim_x.T[:,indx],sim_y.T[:,indy],sim_o.T[:,indo])).A)

    1698. 

  1698.             P.title(str(tlen)+' simulated periods, '+mname)

    1699. 

  1699.             P.xlabel('Time')

    1700. 

  1700.             P.ylabel('Percentage deviation from SS')

    1701. 

  1701.             P.grid()

    1702. 

  1702.             P.legend(leg)

    1703. 

  1703.         elif not indx and indy and indo:

    1704. 

  1704.             for y in indy:

    1705. 

  1705.                 leg.append(conli[y])

    1706. 

  1706.             for o in indo:

    1707. 

  1707.                 leg.append(otherli[o])

    1708. 

  1708.             leg = tuple(leg)

    1709. 

  1709.             figo = P.figure()

    1710. 

  1710.             P.plot(MAT.hstack((sim_y.T[:,indy],sim_o.T[:,indo])).A)

    1711. 

  1711.             P.title(str(tlen)+' simulated periods, '+mname)

    1712. 

  1712.             P.xlabel('Time')

    1713. 

  1713.             P.ylabel('Percentage deviation from SS')

    1714. 

  1714.             P.grid()

    1715. 

  1715.             P.legend(leg)

    1716. 

  1716.         elif indx and not indy and indo:

    1717. 

  1717.             for x in indx:

    1718. 

  1718.                 leg.append(stateli[x])

    1719. 

  1719.             for o in indo:

    1720. 

  1720.                 leg.append(otherli[o])

    1721. 

  1721.             leg = tuple(leg)

    1722. 

  1722.             figo = P.figure()

    1723. 

  1723.             P.plot(MAT.hstack((sim_x.T[:,indx],sim_o.T[:,indo])).A)

    1724. 

  1724.             P.title(str(tlen)+' simulated periods, '+mname)

    1725. 

  1725.             P.xlabel('Time')

    1726. 

  1726.             P.ylabel('Percentage deviation from SS')

    1727. 

  1727.             P.grid()

    1728. 

  1728.             P.legend(leg)

    1729. 

  1729.         elif indx and indy and not indo:

    1730. 

  1730.             for x in indx:

    1731. 

  1731.                 leg.append(stateli[x])

    1732. 

  1732.             for y in indy:

    1733. 

  1733.                 leg.append(conli[y])

    1734. 

  1734.             leg = tuple(leg)

    1735. 

  1735.             figo = P.figure()

    1736. 

  1736.             P.plot(MAT.hstack((sim_x.T[:,indx],sim_y.T[:,indy])).A)

    1737. 

  1737.             P.title(str(tlen)+' simulated periods, '+mname)

    1738. 

  1738.             P.xlabel('Time')

    1739. 

  1739.             P.ylabel('Percentage deviation from SS')

    1740. 

  1740.             P.grid()

    1741. 

  1741.             P.legend(leg)

    1742. 

  1742.         elif indx and not indy and not indo:

    1743. 

  1743.             for x in indx:

    1744. 

  1744.                 leg.append(stateli[x])

    1745. 

  1745.             leg = tuple(leg)

    1746. 

  1746.             figo = P.figure()

    1747. 

  1747.             P.plot(sim_x.T[:,indx].A)

    1748. 

  1748.             P.title(str(tlen)+' simulated periods, '+mname)

    1749. 

  1749.             P.xlabel('Time')

    1750. 

  1750.             P.ylabel('Percentage deviation from SS')

    1751. 

  1751.             P.grid()

    1752. 

  1752.             P.legend(leg)

    1753. 

  1753.         elif not indx and indy and not indo:

    1754. 

  1754.             for y in indy:

    1755. 

  1755.                 leg.append(conli[y])

    1756. 

  1756.             leg = tuple(leg)

    1757. 

  1757.             figo = P.figure()

    1758. 

  1758.             P.plot(sim_y.T[:,indy].A)

    1759. 

  1759.             P.title(str(tlen)+' simulated periods, '+mname)

    1760. 

  1760.             P.xlabel('Time')

    1761. 

  1761.             P.ylabel('Percentage deviation from SS')

    1762. 

  1762.             P.grid()

    1763. 

  1763.             P.legend(leg)

    1764. 

  1764.         elif not indx and not indy and indo:

    1765. 

  1765.             for o in indo:

    1766. 

  1766.                 leg.append(otherli[o])

    1767. 

  1767.             leg = tuple(leg)

    1768. 

  1768.             figo = P.figure()

    1769. 

  1769.             P.plot(sim_o.T[:,indo].A)

    1770. 

  1770.             P.title(str(tlen)+' simulated periods, '+mname)

    1771. 

  1771.             P.xlabel('Time')

    1772. 

  1772.             P.ylabel('Percentage deviation from SS')

    1773. 

  1773.             P.grid()

    1774. 

  1774.             P.legend(leg)

    1775. 

  1775.         return figo

    1776. 

  1776. 

    1777. 

  1777.     def save_sim(self,intup,fpath=None,filtup=None,insim='insim'):

    1778. 

  1778.         # Check if simulations have been carried out

    1779. 

  1779.         if insim not in dir(self):

    1780. 

  1780.             return 'Error: You have not produced any simulations yet! Nothing to graph!'

    1781. 

  1781.         if fpath == None:

    1782. 

  1782.             return 'Error: You have not supplied the function with an fpath= argument where to save the plot!'

    1783. 

  1783.         mname = self.modname

    1784. 

  1784.         vardic = self.vardic

    1785. 

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

    1786. 

  1786.         sim_x = COP.deepcopy(insim[0])

    1787. 

  1787.         sim_y = COP.deepcopy(insim[1])

    1788. 

  1788.         tlen = sim_x.shape[1]

    1789. 

  1789.         if self.oswitch:

    1790. 

  1790.             sim_o = COP.deepcopy(insim[2])

    1791. 

  1791.             nother = self.nother

    1792. 

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

    1793. 

  1793.         nexo = self.nexo

    1794. 

  1794.         nendo = self.nendo

    1795. 

  1795.         ncon = self.ncon

    1796. 

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

    1797. 

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

    1798. 

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

    1799. 

  1799. 

    1800. 

  1800.         # If filtup is None, create filtup from vardic

    1801. 

  1801.         if filtup == None:

    1802. 

  1802.             filli = [0]*len(intup)

    1803. 

  1803.             for i1,x in enumerate(intup):

    1804. 

  1804.                 for y in vardic.keys():

    1805. 

  1805.                     if x in [z[1] for z in vardic[y]['var']]:

    1806. 

  1806.                         indx = [z[1] for z in vardic[y]['var']].index(x)

    1807. 

  1807.                         if 'hp' in [z for z in vardic[y]['mod']][indx]:

    1808. 

  1808.                             filli[i1] = 1

    1809. 

  1809.                         elif 'bk' in [z for z in vardic[y]['mod']][indx]:

    1810. 

  1810.                             filli[i1] = 2

    1811. 

  1811.             filtup = tuple(filli)

    1812. 

  1812. 

    1813. 

  1813. 

    1814. 

  1814.         stateli = exoli+endoli

    1815. 

  1815.         alli = stateli+conli

    1816. 

  1816.         if self.oswitch:

    1817. 

  1817.             alli = alli+otherli

    1818. 

  1818.         # Check if all name in intup are available to graph

    1819. 

  1819.         for name in intup:

    1820. 

  1820.             if name not in alli:

    1821. 

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

    1822. 

  1822.         # Create x, y and z indeces

    1823. 

  1823.         indx = []

    1824. 

  1824.         indy = []

    1825. 

  1825.         indo = []

    1826. 

  1826.         for name in intup:

    1827. 

  1827.             if name in stateli:

    1828. 

  1828.                 indx.append(stateli.index(name))

    1829. 

  1829.             elif name in conli:

    1830. 

  1830.                 indy.append(conli.index(name))

    1831. 

  1831.             elif self.oswitch and name in otherli:

    1832. 

  1832.                 indo.append(otherli.index(name))

    1833. 

  1833.         leg = []

    1834. 

  1834.         indx.sort()

    1835. 

  1835.         indy.sort()

    1836. 

  1836.         indo.sort()     

    1837. 

  1837. 

    1838. 

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

    1839. 

  1839.         for i1 in xrange(sim_y.shape[0]):

    1840. 

  1840.             if indy and (i1 in indy) and filtup[list(intup).index(conli[i1])]:

    1841. 

  1841.                 if filtup[list(intup).index(conli[i1])] == 1:

    1842. 

  1842.                     conli[i1] = conli[i1]+'(hpf)'

    1843. 

  1843.                     yy = sim_y[i1,:].__array__().T

    1844. 

  1844.                     woy = np.zeros((tlen,3))

    1845. 

  1845.                     lam = 1600

    1846. 

  1846.                     yyf = MAT.matrix(hpfilter(data=yy,lam=1600))

    1847. 

  1847.                     sim_y[i1,:] = yyf[0]*100

    1848. 

  1848.                 elif filtup[list(intup).index(conli[i1])] == 2:

    1849. 

  1849.                     conli[i1] = conli[i1]+'(bkf)'

    1850. 

  1850.                     yy = sim_y[i1,:].__array__().T

    1851. 

  1851.                     woy = np.zeros((tlen,1))

    1852. 

  1852.                     up = 6

    1853. 

  1853.                     dn = 32

    1854. 

  1854.                     kkl = 12

    1855. 

  1855.                     yyf = MAT.matrix(bkfilter(yy,up,dn,kkl,tlen))

    1856. 

  1856.                     sim_y[i1,:] = yyf[0]*100

    1857. 

  1857.         # Now filter the state variables!

    1858. 

  1858.         for i1 in xrange(sim_x.shape[0]):

    1859. 

  1859.             if indx and (i1 in indx) and filtup[list(intup).index(stateli[i1])]:

    1860. 

  1860.                 if filtup[list(intup).index(stateli[i1])] == 1:

    1861. 

  1861.                     stateli[i1] = stateli[i1]+'(hpf)'

    1862. 

  1862.                     xx = sim_x[i1,:].__array__().T

    1863. 

  1863.                     wox = np.zeros((tlen,3))

    1864. 

  1864.                     lam = 1600

    1865. 

  1865.                     xxf = MAT.matrix(hpfilter(data=xx,lam=1600))

    1866. 

  1866.                     sim_x[i1,:] = xxf[0]*100

    1867. 

  1867.                 elif filtup[list(intup).index(stateli[i1])] == 2:

    1868. 

  1868.                     stateli[i1] = stateli[i1]+'(bkf)'

    1869. 

  1869.                     xx = sim_x[i1,:].__array__().T

    1870. 

  1870.                     wox = np.zeros((tlen,1))

    1871. 

  1871.                     up = 6

    1872. 

  1872.                     dn = 32

    1873. 

  1873.                     kkl = 12

    1874. 

  1874.                     xxf = MAT.matrix(bkfilter(xx,up,dn,kkl,tlen))

    1875. 

  1875.                     sim_x[i1,:] = xxf[0]*100                  

    1876. 

  1876.             # Now hp filter the other variables!

    1877. 

  1877.         if self.oswitch:

    1878. 

  1878.             for i1 in xrange(sim_o.shape[0]):

    1879. 

  1879.                 if indo and (i1 in indo) and filtup[list(intup).index(otherli[i1])]:

    1880. 

  1880.                     if filtup[list(intup).index(otherli[i1])] == 1:                   

    1881. 

  1881.                         otherli[i1] = otherli[i1]+'(hpf)'

    1882. 

  1882.                         oo = sim_o[i1,:].__array__().T

    1883. 

  1883.                         woo = np.zeros((tlen,3))

    1884. 

  1884.                         lam = 1600

    1885. 

  1885.                         oof = MAT.matrix(hpfilter(data=oo,lam=1600))

    1886. 

  1886.                     elif filtup[list(intup).index(otherli[i1])] == 2:

    1887. 

  1887.                         otherli[i1] = otherli[i1]+'(bkf)'

    1888. 

  1888.                         oo = sim_o[i1,:].__array__().T

    1889. 

  1889.                         woo = np.zeros((tlen,1))

    1890. 

  1890.                         up = 6

    1891. 

  1891.                         dn = 32

    1892. 

  1892.                         kkl = 12

    1893. 

  1893.                         oof = MAT.matrix(bkfilter(oo,up,dn,kkl,tlen)) 

    1894. 

  1894. 

    1895. 

  1895.         if indx and indy and indo:

    1896. 

  1896.             for x in indx:

    1897. 

  1897.                 leg.append(stateli[x])

    1898. 

  1898.             for y in indy:

    1899. 

  1899.                 leg.append(conli[y])

    1900. 

  1900.             for o in indo:

    1901. 

  1901.                 leg.append(otherli[o])

    1902. 

  1902.             leg = tuple(leg)

    1903. 

  1903.             figo = P.figure()

    1904. 

  1904.             P.plot(MAT.hstack((sim_x.T[:,indx],sim_y.T[:,indy],sim_o.T[:,indo])).A)

    1905. 

  1905.             P.title(str(tlen)+' simulated periods, '+mname)

    1906. 

  1906.             P.xlabel('Time')

    1907. 

  1907.             P.ylabel('Percentage deviation from SS')

    1908. 

  1908.             P.grid()

    1909. 

  1909.             P.legend(leg)

    1910. 

  1910.         elif not indx and indy and indo:

    1911. 

  1911.             for y in indy:

    1912. 

  1912.                 leg.append(conli[y])

    1913. 

  1913.             for o in indo:

    1914. 

  1914.                 leg.append(otherli[o])

    1915. 

  1915.             leg = tuple(leg)

    1916. 

  1916.             figo = P.figure()

    1917. 

  1917.             P.plot(MAT.hstack((sim_y.T[:,indy],sim_o.T[:,indo])).A)

    1918. 

  1918.             P.title(str(tlen)+' simulated periods, '+mname)

    1919. 

  1919.             P.xlabel('Time')

    1920. 

  1920.             P.ylabel('Percentage deviation from SS')

    1921. 

  1921.             P.grid()

    1922. 

  1922.             P.legend(leg)

    1923. 

  1923.         elif indx and not indy and indo:

    1924. 

  1924.             for x in indx:

    1925. 

  1925.                 leg.append(stateli[x])

    1926. 

  1926.             for o in indo:

    1927. 

  1927.                 leg.append(otherli[o])

    1928. 

  1928.             leg = tuple(leg)

    1929. 

  1929.             figo = P.figure()

    1930. 

  1930.             P.plot(MAT.hstack((sim_x.T[:,indx],sim_o.T[:,indo])).A)

    1931. 

  1931.             P.title(str(tlen)+' simulated periods, '+mname)

    1932. 

  1932.             P.xlabel('Time')

    1933. 

  1933.             P.ylabel('Percentage deviation from SS')

    1934. 

  1934.             P.grid()

    1935. 

  1935.             P.legend(leg)

    1936. 

  1936.         elif indx and indy and not indo:

    1937. 

  1937.             for x in indx:

    1938. 

  1938.                 leg.append(stateli[x])

    1939. 

  1939.             for y in indy:

    1940. 

  1940.                 leg.append(conli[y])

    1941. 

  1941.             leg = tuple(leg)

    1942. 

  1942.             figo = P.figure()

    1943. 

  1943.             P.plot(MAT.hstack((sim_x.T[:,indx],sim_y.T[:,indy])).A)

    1944. 

  1944.             P.title(str(tlen)+' simulated periods, '+mname)

    1945. 

  1945.             P.xlabel('Time')

    1946. 

  1946.             P.ylabel('Percentage deviation from SS')

    1947. 

  1947.             P.grid()

    1948. 

  1948.             P.legend(leg)

    1949. 

  1949.         elif indx and not indy and not indo:

    1950. 

  1950.             for x in indx:

    1951. 

  1951.                 leg.append(stateli[x])

    1952. 

  1952.             leg = tuple(leg)

    1953. 

  1953.             figo = P.figure()

    1954. 

  1954.             P.plot(sim_x.T[:,indx].A)

    1955. 

  1955.             P.title(str(tlen)+' simulated periods, '+mname)

    1956. 

  1956.             P.xlabel('Time')

    1957. 

  1957.             P.ylabel('Percentage deviation from SS')

    1958. 

  1958.             P.grid()

    1959. 

  1959.             P.legend(leg)

    1960. 

  1960.         elif not indx and indy and not indo:

    1961. 

  1961.             for y in indy:

    1962. 

  1962.                 leg.append(conli[y])

    1963. 

  1963.             leg = tuple(leg)

    1964. 

  1964.             figo = P.figure()

    1965. 

  1965.             P.plot(sim_y.T[:,indy].A)

    1966. 

  1966.             P.title(str(tlen)+' simulated periods, '+mname)

    1967. 

  1967.             P.xlabel('Time')

    1968. 

  1968.             P.ylabel('Percentage deviation from SS')

    1969. 

  1969.             P.grid()

    1970. 

  1970.             P.legend(leg)

    1971. 

  1971.         elif not indx and not indy and indo:

    1972. 

  1972.             for o in indo:

    1973. 

  1973.                 leg.append(otherli[o])

    1974. 

  1974.             leg = tuple(leg)

    1975. 

  1975.             figo = P.figure()

    1976. 

  1976.             P.plot(sim_o.T[:,indo].A)

    1977. 

  1977.             P.title(str(tlen)+' simulated periods, '+mname)

    1978. 

  1978.             P.xlabel('Time')

    1979. 

  1979.             P.ylabel('Percentage deviation from SS')

    1980. 

  1980.             P.grid()

    1981. 

  1981.             P.legend(leg)

    1982. 

  1982.         PLT.savefig(fpath,dpi=100)

    1983. 

  1983.         return "Your plot has been saved in: "+fpath

    1984. 

  1984. 

    1985. 

  1985. 

    1986. 

  1986.     def irf(self,tlen,sntup,shockmod=None):

    1987. 

  1987.         # Deals with 0-tuples

    1988. 

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

    1989. 

  1989.             sntup = (sntup,)

    1990. 

  1990.         # Deal with shockmod, also take care of 0-tuples by turning all into lists

    1991. 

  1991.         if shockmod != None and type(shockmod) != type((1,2,3)) and type(shockmod) == type(1.0):

    1992. 

  1992.             shockmod = [shockmod,]

    1993. 

  1993.         elif shockmod != None and len(shockmod) > 1:

    1994. 

  1994.             shockmod = [x for x in shockmod]

    1995. 

  1995.         elif shockmod == None:

    1996. 

  1996.             shockmod = [1.0,]*len(sntup)

    1997. 

  1997.         tlen = tlen + 1

    1998. 

  1998.         ncon = self.ncon

    1999. 

  1999.         nexo = self.nexo

    2000. 

  2000.         nendo = self.nendo

    2001. 

  2001.         tstates = self.tstates

    2002. 

  2002.         sigma = self.sigma

    2003. 

  2003.         ssigma = self.ssigma

    2004. 

  2004.         if self.oswitch:

    2005. 

  2005.             numjl = self.numjl

    2006. 

  2006.             numhl = self.numhl

    2007. 

  2007.         kx = self.KX

    2008. 

  2008.         ky = self.KY

    2009. 

  2009.         pp = self.P

    2010. 

  2010.         ff = self.F

    2011. 

  2011.         ee = self.E

    2012. 

  2012.         gg = self.G

    2013. 

  2013. 

    2014. 

  2014.         sposli=[]

    2015. 

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

    2016. 

  2016.         # Check if names are valid

    2017. 

  2017.         for name in sntup:

    2018. 

  2018.             if name not in exoli:

    2019. 

  2019.                 return 'Error: '+name+' is not a valid exoshock name for this model!'

    2020. 

  2020.         for name in sntup:

    2021. 

  2021.             sposli.append(exoli.index(name))

    2022. 

  2022. 

    2023. 

  2023.         # Expose spos choice to self for show_irf

    2024. 

  2024.         self.spos = (sntup,sposli)

    2025. 

  2025. 

    2026. 

  2026.         shock = MAT.zeros((nexo,1))

    2027. 

  2027.         sendo = MAT.zeros((nendo,1))

    2028. 

  2028. 

    2029. 

  2029.         for i1,spos in enumerate(sposli):

    2030. 

  2030.             shock[spos,0] = 1.0*shockmod[i1]

    2031. 

  2031.         shock = MAT.vstack((shock,sendo))

    2032. 

  2032.         

    2033. 

  2033.         self.sshock = COP.deepcopy(shock)

    2034. 

  2034. 

    2035. 

  2035.         x_one_m1 = shock

    2036. 

  2036.         x_one_0 = pp*x_one_m1

    2037. 

  2037.         x_two_m1 = shock

    2038. 

  2038.         y_one_0 = ff*x_one_m1

    2039. 

  2039.         # Because they are jump variables

    2040. 

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

    2041. 

  2041.         y_two_0 = 0.5*MAT.kron(MAT.eye(ncon),x_one_m1.T)*ee*x_one_m1

    2042. 

  2042.         if self.oswitch:

    2043. 

  2043.             nother = self.nother

    2044. 

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

    2045. 

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

    2046. 

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

    2047. 

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

    2048. 

  2048. 

    2049. 

  2049.         x_one_c = COP.deepcopy(x_one_m1)

    2050. 

  2050.         y_one_c = COP.deepcopy(y_one_0)

    2051. 

  2051.         x_two_c = COP.deepcopy(x_two_m1)

    2052. 

  2052.         y_two_c = COP.deepcopy(y_two_0)

    2053. 

  2053.         if self.oswitch:

    2054. 

  2054.             o_one_c = COP.deepcopy(o_one_0)

    2055. 

  2055.             o_two_c = COP.deepcopy(o_two_0)

    2056. 

  2056. 

    2057. 

  2057.         x_one = COP.deepcopy(x_one_c)

    2058. 

  2058.         y_one = COP.deepcopy(y_one_c)

    2059. 

  2059.         x_two = COP.deepcopy(x_two_c)

    2060. 

  2060.         y_two = COP.deepcopy(y_two_c)

    2061. 

  2061.         if self.oswitch:

    2062. 

  2062.             o_one = COP.deepcopy(o_one_c)

    2063. 

  2063.             o_two = COP.deepcopy(o_two_c)

    2064. 

  2064. 

    2065. 

  2065.         for i1 in range(tlen):

    2066. 

  2066.             x_one_n = pp*x_one_c

    2067. 

  2067.             y_one_n = ff*x_one_c

    2068. 

  2068.             x_two_n = pp*x_two_c+0.5*MAT.kron(MAT.eye(tstates),x_one_c.T)*gg*x_one_c

    2069. 

  2069.             y_two_n = ff*x_two_c+0.5*MAT.kron(MAT.eye(ncon),x_one_c.T)*ee*x_one_c

    2070. 

  2070.             if self.oswitch:

    2071. 

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

    2072. 

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

    2073. 

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

    2074. 

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

    2075. 

  2075. 

    2076. 

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

    2077. 

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

    2078. 

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

    2079. 

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

    2080. 

  2080.             if self.oswitch:

    2081. 

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

    2082. 

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

    2083. 

  2083. 

    2084. 

  2084.             x_one_c = COP.deepcopy(x_one_n)

    2085. 

  2085.             y_one_c = COP.deepcopy(y_one_n)

    2086. 

  2086.             x_two_c = COP.deepcopy(x_one_n)

    2087. 

  2087.             y_two_c = COP.deepcopy(y_two_n)

    2088. 

  2088. 

    2089. 

  2089.         self.irf_x_two = x_two[:,1:-1]

    2090. 

  2090.         self.irf_y_two = y_two[:,1:-1]

    2091. 

  2091.         self.inirf = [self.irf_x_two,self.irf_y_two]

    2092. 

  2092.         if self.oswitch:

    2093. 

  2093.             self.irf_o_one = o_one[:,2:]

    2094. 

  2094.             self.irf_o_two = o_two[:,2:]

    2095. 

  2095.             self.inirf = self.inirf + [self.irf_o_two,]

    2096. 

  2096.         

    2097. 

  2097. 

    2098. 

  2098.     def show_irf(self,intup,inirf='inirf'):

    2099. 

  2099.         # Check if simulations have been carried out

    2100. 

  2100.         if inirf not in dir(self):

    2101. 

  2101.             return 'Error: You have not produced any IRFs yet! Nothing to graph!'

    2102. 

  2102.         inirf = eval('self.'+inirf)

    2103. 

  2103.         irf_x = COP.deepcopy(inirf[0])

    2104. 

  2104.         tlen = irf_x.shape[1]

    2105. 

  2105.         irf_x = MAT.hstack((MAT.zeros((irf_x.shape[0],20)),irf_x))

    2106. 

  2106.         irf_y = COP.deepcopy(inirf[1])

    2107. 

  2107.         irf_y = MAT.hstack((MAT.zeros((irf_y.shape[0],20)),irf_y))

    2108. 

  2108.         #irf_x = irf_x*100.0

    2109. 

  2109.         #irf_y = irf_y*100.0

    2110. 

  2110.         mname = self.modname

    2111. 

  2111.         vardic = self.vardic

    2112. 

  2112.         time_axis = np.arange(-20,tlen,1)

    2113. 

  2113.         if self.oswitch:

    2114. 

  2114.             irf_o = COP.deepcopy(inirf[2])

    2115. 

  2115.             irf_o = MAT.hstack((MAT.zeros((irf_o.shape[0],20)),irf_o))

    2116. 

  2116.             nother = self.nother

    2117. 

  2117.             varother = self.vardic['other']['var']

    2118. 

  2118.             otherli = [x[1] for x in varother]

    2119. 

  2119.         nexo = self.nexo

    2120. 

  2120.         nendo = self.nendo

    2121. 

  2121.         ncon = self.ncon

    2122. 

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

    2123. 

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

    2124. 

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

    2125. 

  2125. 

    2126. 

  2126.         stateli = exoli+endoli

    2127. 

  2127.         alli = stateli+conli

    2128. 

  2128.         if self.oswitch:

    2129. 

  2129.             alli = alli+otherli

    2130. 

  2130.         # Check if all name in intup are available to graph

    2131. 

  2131.         for name in intup:

    2132. 

  2132.             if name not in alli:

    2133. 

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

    2134. 

  2134.         # Create x, y and z indeces

    2135. 

  2135.         indx = []

    2136. 

  2136.         indy = []

    2137. 

  2137.         indo = []

    2138. 

  2138.         for name in intup:

    2139. 

  2139.             if name in stateli:

    2140. 

  2140.                 indx.append(stateli.index(name))

    2141. 

  2141.             elif name in conli:

    2142. 

  2142.                 indy.append(conli.index(name))

    2143. 

  2143.             elif self.oswitch and name in otherli:

    2144. 

  2144.                 indo.append(otherli.index(name))

    2145. 

  2145.         leg = []

    2146. 

  2146.         indx.sort()

    2147. 

  2147.         indy.sort()

    2148. 

  2148.         indo.sort()

    2149. 

  2149. 

    2150. 

  2150.         if indx and indy and indo:

    2151. 

  2151.             for x in indx:

    2152. 

  2152.                 leg.append(stateli[x])

    2153. 

  2153.             for y in indy:

    2154. 

  2154.                 leg.append(conli[y])

    2155. 

  2155.             for o in indo:

    2156. 

  2156.                 leg.append(otherli[o])

    2157. 

  2157.             leg = tuple(leg)

    2158. 

  2158.             figo = P.figure()

    2159. 

  2159.             ax = figo.add_subplot(111)

    2160. 

  2160.             ax.set_xlim([-20,tlen])

    2161. 

  2161.             ax.plot(time_axis,MAT.hstack((irf_x.T[:,indx],irf_y.T[:,indy],irf_o.T[:,indo])).A)

    2162. 

  2162.             P.title(str(tlen)+' IRF periods, '+mname)

    2163. 

  2163.             P.xlabel('Time')

    2164. 

  2164.             P.ylabel('Percentage Deviation from SS')

    2165. 

  2165.             P.grid()

    2166. 

  2166.             P.legend(leg)

    2167. 

  2167.         elif not indx and indy and indo:

    2168. 

  2168.             for y in indy:

    2169. 

  2169.                 leg.append(conli[y])

    2170. 

  2170.             for o in indo:

    2171. 

  2171.                 leg.append(otherli[o])

    2172. 

  2172.             leg = tuple(leg)

    2173. 

  2173.             figo = P.figure()

    2174. 

  2174.             ax = figo.add_subplot(111)

    2175. 

  2175.             ax.set_xlim([-20,tlen])

    2176. 

  2176.             ax.plot(time_axis,MAT.hstack((irf_y.T[:,indy],irf_o.T[:,indo])).A)

    2177. 

  2177.             P.title(str(tlen)+' IRF periods, '+mname)

    2178. 

  2178.             P.xlabel('Time')

    2179. 

  2179.             P.ylabel('Percentage Deviation from SS')

    2180. 

  2180.             P.grid()

    2181. 

  2181.             P.legend(leg)

    2182. 

  2182.         elif indx and not indy and indo:

    2183. 

  2183.             for x in indx:

    2184. 

  2184.                 leg.append(stateli[x])

    2185. 

  2185.             for o in indo:

    2186. 

  2186.                 leg.append(otherli[o])

    2187. 

  2187.             leg = tuple(leg)

    2188. 

  2188.             figo = P.figure()

    2189. 

  2189.             ax = figo.add_subplot(111)

    2190. 

  2190.             ax.set_xlim([-20,tlen])

    2191. 

  2191.             ax.plot(time_axis,MAT.hstack((irf_x.T[:,indx],irf_o.T[:,indo])).A)

    2192. 

  2192.             P.title(str(tlen)+' IRF periods, '+mname)

    2193. 

  2193.             P.xlabel('Time')

    2194. 

  2194.             P.ylabel('Percentage Deviation from SS')

    2195. 

  2195.             P.grid()

    2196. 

  2196.             P.legend(leg)

    2197. 

  2197.         elif indx and indy and not indo:

    2198. 

  2198.             for x in indx:

    2199. 

  2199.                 leg.append(stateli[x])

    2200. 

  2200.             for y in indy:

    2201. 

  2201.                 leg.append(conli[y])

    2202. 

  2202.             leg = tuple(leg)

    2203. 

  2203.             figo = P.figure()

    2204. 

  2204.             ax = figo.add_subplot(111)

    2205. 

  2205.             ax.set_xlim([-20,tlen])

    2206. 

  2206.             ax.plot(time_axis,MAT.hstack((irf_x.T[:,indx],irf_y.T[:,indy])).A)

    2207. 

  2207.             P.title(str(tlen)+' IRF periods, '+mname)

    2208. 

  2208.             P.xlabel('Time')

    2209. 

  2209.             P.ylabel('Percentage Deviation from SS')

    2210. 

  2210.             P.grid()

    2211. 

  2211.             P.legend(leg)

    2212. 

  2212.         elif indx and not indy and not indo:

    2213. 

  2213.             for x in indx:

    2214. 

  2214.                 leg.append(stateli[x])

    2215. 

  2215.             leg = tuple(leg)

    2216. 

  2216.             figo = P.figure()

    2217. 

  2217.             ax = figo.add_subplot(111)

    2218. 

  2218.             ax.set_xlim([-20,tlen])

    2219. 

  2219.             ax.plot(time_axis,irf_x.T[:,indx].A)

    2220. 

  2220.             P.title(str(tlen)+' IRF periods, '+mname)

    2221. 

  2221.             P.xlabel('Time')

    2222. 

  2222.             P.ylabel('Percentage Deviation from SS')

    2223. 

  2223.             P.grid()

    2224. 

  2224.             P.legend(leg)

    2225. 

  2225.         elif not indx and indy and not indo:

    2226. 

  2226.             for y in indy:

    2227. 

  2227.                 leg.append(conli[y])

    2228. 

  2228.             leg = tuple(leg)

    2229. 

  2229.             figo = P.figure()

    2230. 

  2230.             ax = figo.add_subplot(111)

    2231. 

  2231.             ax.set_xlim([-20,tlen])

    2232. 

  2232.             ax.plot(time_axis,irf_y.T[:,indy].A)

    2233. 

  2233.             P.title(str(tlen)+' IRF periods, '+mname)

    2234. 

  2234.             P.xlabel('Time')

    2235. 

  2235.             P.ylabel('Percentage Deviation from SS')

    2236. 

  2236.             P.grid()

    2237. 

  2237.             P.legend(leg)

    2238. 

  2238.         elif not indx and not indy and indo:

    2239. 

  2239.             for o in indo:

    2240. 

  2240.                 leg.append(otherli[o])

    2241. 

  2241.             leg = tuple(leg)

    2242. 

  2242.             figo = P.figure()

    2243. 

  2243.             ax = figo.add_subplot(111)

    2244. 

  2244.             ax.set_xlim([-20,tlen])

    2245. 

  2245.             ax.plot(time_axis,irf_o.T[:,indo].A)

    2246. 

  2246.             P.title(str(tlen)+' IRF periods, '+mname)

    2247. 

  2247.             P.xlabel('Time')

    2248. 

  2248.             P.ylabel('Percentage Deviation from SS')

    2249. 

  2249.             P.grid()

    2250. 

  2250.             P.legend(leg)

    2251. 

  2251.         return figo

    2252. 

  2252. 

    2253. 

  2253.     def save_irf(self,intup,fpath=None,inirf='inirf'):

    2254. 

  2254.         # Check if simulations have been carried out

    2255. 

  2255.         if inirf not in dir(self):

    2256. 

  2256.             return 'Error: You have not produced any IRFs yet! Nothing to graph!'

    2257. 

  2257.         if fpath == None:

    2258. 

  2258.             return 'Error: You have not supplied the function with an fpath= argument where to save the plot!'

    2259. 

  2259.         inirf = eval('self.'+inirf)

    2260. 

  2260.         irf_x = COP.deepcopy(inirf[0])

    2261. 

  2261.         tlen = irf_x.shape[1]

    2262. 

  2262.         irf_x = MAT.hstack((MAT.zeros((irf_x.shape[0],20)),irf_x))

    2263. 

  2263.         irf_y = COP.deepcopy(inirf[1])

    2264. 

  2264.         irf_y = MAT.hstack((MAT.zeros((irf_y.shape[0],20)),irf_y))

    2265. 

  2265.         irf_x = irf_x*100

    2266. 

  2266.         irf_y = irf_y*100

    2267. 

  2267.         mname = self.modname

    2268. 

  2268.         vardic = self.vardic

    2269. 

  2269.         time_axis = np.arange(-20,tlen,1)

    2270. 

  2270.         if self.oswitch:

    2271. 

  2271.             irf_o = COP.deepcopy(inirf[2])

    2272. 

  2272.             nother = self.nother

    2273. 

  2273.             varother = self.vardic['other']['var']

    2274. 

  2274.             otherli = [x[1] for x in varother]

    2275. 

  2275.         nexo = self.nexo

    2276. 

  2276.         nendo = self.nendo

    2277. 

  2277.         ncon = self.ncon

    2278. 

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

    2279. 

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

    2280. 

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

    2281. 

  2281. 

    2282. 

  2282.         stateli = exoli+endoli

    2283. 

  2283.         alli = stateli+conli

    2284. 

  2284.         if self.oswitch:

    2285. 

  2285.             alli = alli+otherli

    2286. 

  2286.         # Check if all name in intup are available to graph

    2287. 

  2287.         for name in intup:

    2288. 

  2288.             if name not in alli:

    2289. 

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

    2290. 

  2290.         # Create x, y and z indeces

    2291. 

  2291.         indx = []

    2292. 

  2292.         indy = []

    2293. 

  2293.         indo = []

    2294. 

  2294.         for name in intup:

    2295. 

  2295.             if name in stateli:

    2296. 

  2296.                 indx.append(stateli.index(name))

    2297. 

  2297.             elif name in conli:

    2298. 

  2298.                 indy.append(conli.index(name))

    2299. 

  2299.             elif self.oswitch and name in otherli:

    2300. 

  2300.                 indo.append(otherli.index(name))

    2301. 

  2301.         leg = []

    2302. 

  2302.         indx.sort()

    2303. 

  2303.         indy.sort()

    2304. 

  2304.         indo.sort()

    2305. 

  2305. 

    2306. 

  2306.         if indx and indy and indo:

    2307. 

  2307.             for x in indx:

    2308. 

  2308.                 leg.append(stateli[x])

    2309. 

  2309.             for y in indy:

    2310. 

  2310.                 leg.append(conli[y])

    2311. 

  2311.             for o in indo:

    2312. 

  2312.                 leg.append(otherli[o])

    2313. 

  2313.             leg = tuple(leg)

    2314. 

  2314.             figo = P.figure()

    2315. 

  2315.             ax = figo.add_subplot(111)

    2316. 

  2316.             ax.set_xlim([-20,tlen])

    2317. 

  2317.             ax.plot(time_axis,MAT.hstack((irf_x.T[:,indx],irf_y.T[:,indy],irf_o.T[:,indo])).A)

    2318. 

  2318.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2319. 

  2319.             P.xlabel('Time')

    2320. 

  2320.             P.ylabel('Percentage Deviation from SS')

    2321. 

  2321.             P.grid()

    2322. 

  2322.             P.legend(leg)

    2323. 

  2323.         elif not indx and indy and indo:

    2324. 

  2324.             for y in indy:

    2325. 

  2325.                 leg.append(conli[y])

    2326. 

  2326.             for o in indo:

    2327. 

  2327.                 leg.append(otherli[o])

    2328. 

  2328.             leg = tuple(leg)

    2329. 

  2329.             figo = P.figure()

    2330. 

  2330.             ax = figo.add_subplot(111)

    2331. 

  2331.             ax.set_xlim([-20,tlen])

    2332. 

  2332.             ax.plot(time_axis,MAT.hstack((irf_y.T[:,indy],irf_o.T[:,indo])).A)

    2333. 

  2333.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2334. 

  2334.             P.xlabel('Time')

    2335. 

  2335.             P.ylabel('Percentage Deviation from SS')

    2336. 

  2336.             P.grid()

    2337. 

  2337.             P.legend(leg)

    2338. 

  2338.         elif indx and not indy and indo:

    2339. 

  2339.             for x in indx:

    2340. 

  2340.                 leg.append(stateli[x])

    2341. 

  2341.             for o in indo:

    2342. 

  2342.                 leg.append(otherli[o])

    2343. 

  2343.             leg = tuple(leg)

    2344. 

  2344.             figo = P.figure()

    2345. 

  2345.             ax = figo.add_subplot(111)

    2346. 

  2346.             ax.set_xlim([-20,tlen])

    2347. 

  2347.             ax.plot(time_axis,MAT.hstack((irf_x.T[:,indx],irf_o.T[:,indo])).A)

    2348. 

  2348.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2349. 

  2349.             P.xlabel('Time')

    2350. 

  2350.             P.ylabel('Percentage Deviation from SS')

    2351. 

  2351.             P.grid()

    2352. 

  2352.             P.legend(leg)

    2353. 

  2353.         elif indx and indy and not indo:

    2354. 

  2354.             for x in indx:

    2355. 

  2355.                 leg.append(stateli[x])

    2356. 

  2356.             for y in indy:

    2357. 

  2357.                 leg.append(conli[y])

    2358. 

  2358.             leg = tuple(leg)

    2359. 

  2359.             figo = P.figure()

    2360. 

  2360.             ax = figo.add_subplot(111)

    2361. 

  2361.             ax.set_xlim([-20,tlen])

    2362. 

  2362.             ax.plot(time_axis,MAT.hstack((irf_x.T[:,indx],irf_y.T[:,indy])).A)

    2363. 

  2363.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2364. 

  2364.             P.xlabel('Time')

    2365. 

  2365.             P.ylabel('Percentage Deviation from SS')

    2366. 

  2366.             P.grid()

    2367. 

  2367.             P.legend(leg)

    2368. 

  2368.         elif indx and not indy and not indo:

    2369. 

  2369.             for x in indx:

    2370. 

  2370.                 leg.append(stateli[x])

    2371. 

  2371.             leg = tuple(leg)

    2372. 

  2372.             figo = P.figure()

    2373. 

  2373.             ax = figo.add_subplot(111)

    2374. 

  2374.             ax.set_xlim([-20,tlen])

    2375. 

  2375.             ax.plot(time_axis,irf_x.T[:,indx].A)

    2376. 

  2376.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2377. 

  2377.             P.xlabel('Time')

    2378. 

  2378.             P.ylabel('Percentage Deviation from SS')

    2379. 

  2379.             P.grid()

    2380. 

  2380.             P.legend(leg)

    2381. 

  2381.         elif not indx and indy and not indo:

    2382. 

  2382.             for y in indy:

    2383. 

  2383.                 leg.append(conli[y])

    2384. 

  2384.             leg = tuple(leg)

    2385. 

  2385.             figo = P.figure()

    2386. 

  2386.             ax = figo.add_subplot(111)

    2387. 

  2387.             ax.set_xlim([-20,tlen])

    2388. 

  2388.             ax.plot(time_axis,irf_y.T[:,indy].A)

    2389. 

  2389.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2390. 

  2390.             P.xlabel('Time')

    2391. 

  2391.             P.ylabel('Percentage Deviation from SS, hp-filtered')

    2392. 

  2392.             P.grid()

    2393. 

  2393.             P.legend(leg)

    2394. 

  2394.         elif not indx and not indy and indo:

    2395. 

  2395.             for o in indo:

    2396. 

  2396.                 leg.append(otherli[o])

    2397. 

  2397.             leg = tuple(leg)

    2398. 

  2398.             figo = P.figure()

    2399. 

  2399.             ax = figo.add_subplot(111)

    2400. 

  2400.             ax.set_xlim([-20,tlen])

    2401. 

  2401.             ax.plot(time_axis,irf_o.T[:,indo].A)

    2402. 

  2402.             P.title(str(tlen)+' simulated IRF periods, '+mname)

    2403. 

  2403.             P.xlabel('Time')

    2404. 

  2404.             P.ylabel('Percentage Deviation from SS, hp-filtered')

    2405. 

  2405.             P.grid()

    2406. 

  2406.             P.legend(leg)

    2407. 

  2407.         PLT.savefig(fpath,dpi=100)

    2408. 

  2408.         return "Your plot has been save in: "+fpath

    2409. 

  2409. 

    2410. 

  2410. #----------------------------------------------------------------------------------------------------------------------

    2411. 

  2411. class MatKlein2D(PyKlein2D):

    2412. 

  2412. 

    2413. 

  2413.     def __init__(self,intup):

    2414. 

  2414.         self.sess1 = intup[-1]

    2415. 

  2415.         intup2 = intup[:-1]

    2416. 

  2416.         # Get all PyKlein2D attributes

    2417. 

  2417.         PyKlein2D.__init__(self, intup2)

    2418. 

  2418. 

    2419. 

  2419.     def solve(self):

    2420. 

  2420.         tstates = self.tstates

    2421. 

  2421.         sess1 = self.sess1

    2422. 

  2422.         ssigma = self.ssigma

    2423. 

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

    2424. 

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

    2425. 

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

    2426. 

  2426.         mlabraw.put(sess1,'jac',self.gra)

    2427. 

  2427.         mlabraw.put(sess1,'hess',self.hes)

    2428. 

  2428.         mlabraw.put(sess1,'nstates',self.tstates)

    2429. 

  2429.         mlabraw.put(sess1,'ssigma',self.ssigma)

    2430. 

  2430.         mlabraw.eval(sess1,'[ff,pp,ee,gg,kx,ky] = solab2(jac,hess,ssigma,nstates)')

    2431. 

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

    2432. 

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

    2433. 

  2433.         self.E = np.matrix(mlabraw.get(sess1,'ee'))

    2434. 

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

    2435. 

  2435.         self.KX = np.matrix(mlabraw.get(sess1,'kx'))

    2436. 

  2436.         self.KY = np.matrix(mlabraw.get(sess1,'ky'))

    2437. 

  2437. #----------------------------------------------------------------------------------------------------------------------

    2438. 

  2438. class ForKleinD(PyKlein2D):

    2439. 

  2439. 

    2440. 

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

    2441. 

  2441.         if other != None:

    2442. 

  2442.             self.other = other

    2443. 

  2443.         self.gra = intup[0]

    2444. 

  2444.         self.nendo = intup[1]

    2445. 

  2445.         self.nexo = intup[2]

    2446. 

  2446.         self.ncon = intup[3]

    2447. 

  2447.         self.sigma = intup[4]

    2448. 

  2448.         self.A = intup[5]

    2449. 

  2449.         self.B = intup[6]

    2450. 

  2450.         self.vardic = intup[7]

    2451. 

  2451.         self.vdic = intup[8]

    2452. 

  2452.         self.modname = intup[9]

    2453. 

  2453.         self.audic = intup[10]

    2454. 

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

    2455. 

  2455.             self.numjl = intup[11]

    2456. 

  2456.             self.nother = intup[12]

    2457. 

  2457.             self.oswitch = 1

    2458. 

  2458.         else:

    2459. 

  2459.             self.oswitch = 0

    2460. 

  2460.         self.tstates = self.nendo+self.nexo

    2461. 

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

    2462. 

  2462. 

    2463. 

  2463.     def solve(self):

    2464. 

  2464.         A = self.A

    2465. 

  2465.         B = self.B

    2466. 

  2466.         tstates = self.tstates

    2467. 

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

    2468. 

  2468.         self.alpha = abs(alpha)

    2469. 

  2469.         self.betta = abs(betta)

    2470. 

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

    2471. 

  2471.             return

    2472. 

  2472.         else:

    2473. 

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

    2474. 

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

    2475. 

  2475. 

    2476. 

  2476. 

    2477. 

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

    2478. 

  2478.         # Deals with 0-tuples

    2479. 

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

    2480. 

  2480.             sntup = (sntup,)

    2481. 

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

    2482. 

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

    2483. 

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

    2484. 

  2484.         indx=[]

    2485. 

  2485.         if sntup == None:

    2486. 

  2486.             indx = range(len(exoli))

    2487. 

  2487.         else:

    2488. 

  2488.             for name in sntup:

    2489. 

  2489.                 if name not in exoli:

    2490. 

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

    2491. 

  2491.                 else:

    2492. 

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

    2493. 

  2493.         indx.sort()

    2494. 

  2494.         ncon = self.ncon

    2495. 

  2495.         nexo = self.nexo

    2496. 

  2496.         nendo = self.nendo

    2497. 

  2497.         tstates = self.tstates

    2498. 

  2498.         tlena = 1000+tlen

    2499. 

  2499.         sigma = self.sigma

    2500. 

  2500.         ssigma = self.ssigma

    2501. 

  2501.         if self.oswitch:

    2502. 

  2502.             numjl = self.numjl

    2503. 

  2503. 

    2504. 

  2504.         pp = self.P

    2505. 

  2505.         ff = self.F

    2506. 

  2506. 

    2507. 

  2507.         if shockvec == None:

    2508. 

  2508.             count = 0

    2509. 

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

    2510. 

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

    2511. 

  2511.                     if count in indx:

    2512. 

  2512.                         ranvec = MAT.vstack((ranvec,np.sqrt(varia)*MAT.matrix(np.random.standard_normal(tlena))))

    2513. 

  2513.                     else:

    2514. 

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

    2515. 

  2515.                 else:

    2516. 

  2516.                     if count in indx:

    2517. 

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

    2518. 

  2518.                     else:

    2519. 

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

    2520. 

  2520.                 count = count + 1

    2521. 

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

    2522. 

  2522.             # Save the random shocks vector

    2523. 

  2523.             self.shockvec = COP.deepcopy(ranvec)

    2524. 

  2524.         else:

    2525. 

  2525.             ranvec = COP.deepcopy(shockvec)

    2526. 

  2526. 

    2527. 

  2527.         x_one_m1 = ranvec[:,0]

    2528. 

  2528.         y_one_0 = ff*x_one_m1

    2529. 

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

    2530. 

  2530.         x_one_0 = pp*x_one_m1

    2531. 

  2531.         if self.oswitch:

    2532. 

  2532.             nother = self.nother

    2533. 

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

    2534. 

  2534. 

    2535. 

  2535.         x_one_c = COP.deepcopy(x_one_m1)

    2536. 

  2536.         y_one_c = COP.deepcopy(y_one_0)

    2537. 

  2537.         if self.oswitch:

    2538. 

  2538.             o_one_c = COP.deepcopy(o_one_0)

    2539. 

  2539. 

    2540. 

  2540.         x_one = COP.deepcopy(x_one_c)

    2541. 

  2541.         y_one = COP.deepcopy(y_one_c)

    2542. 

  2542.         if self.oswitch:

    2543. 

  2543.             o_one = COP.deepcopy(o_one_c)

    2544. 

  2544. 

    2545. 

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

    2546. 

  2546. 

    2547. 

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

    2548. 

  2548.             y_one_n = ff*x_one_c

    2549. 

  2549.             if self.oswitch:

    2550. 

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

    2551. 

  2551. 

    2552. 

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

    2553. 

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

    2554. 

  2554.             if self.oswitch:

    2555. 

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

    2556. 

  2556. 

    2557. 

  2557.             x_one_c = x_one_n

    2558. 

  2558.             y_one_c = y_one_n

    2559. 

  2559. 

    2560. 

  2560.         # Throw away first 1000 obs

    2561. 

  2561.         x_one = x_one[:,1000:]

    2562. 

  2562.         y_one = y_one[:,1000:]

    2563. 

  2563.         if self.oswitch:

    2564. 

  2564.             o_one = o_one[:,1000:]

    2565. 

  2565. 

    2566. 

  2566.         self.sim_x_one = x_one

    2567. 

  2567.         self.sim_y_one = y_one

    2568. 

  2568.         self.insim = [self.sim_x_one,self.sim_y_one]

    2569. 

  2569.         if self.oswitch:

    2570. 

  2570.             self.sim_o_one = o_one

    2571. 

  2571.             self.insim = self.insim + [self.sim_o_one,]

    2572. 

  2572. 

    2573. 

  2573.     def irf(self,tlen,sntup=None,shockmod=None):

    2574. 

  2574.         # Deals with 0-tuples

    2575. 

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

    2576. 

  2576.             sntup = (sntup,)

    2577. 

  2577.         # Deal with shockmod, also take care of 0-tuples by turning all into lists

    2578. 

  2578.         if shockmod != None and type(shockmod) != type((1,2,3)) and type(shockmod) == type(1.0):

    2579. 

  2579.             shockmod = [shockmod,]

    2580. 

  2580.         elif shockmod != None and len(shockmod) > 1:

    2581. 

  2581.             shockmod = [x for x in shockmod]

    2582. 

  2582.         elif shockmod == None:

    2583. 

  2583.             shockmod = [1.0,]*len(sntup)

    2584. 

  2584.         tlen = tlen + 1

    2585. 

  2585.         ncon = self.ncon

    2586. 

  2586.         nexo = self.nexo

    2587. 

  2587.         nendo = self.nendo

    2588. 

  2588.         tstates = self.tstates

    2589. 

  2589.         sigma = self.sigma

    2590. 

  2590.         ssigma = self.ssigma

    2591. 

  2591.         tlen = tlen

    2592. 

  2592.         if self.oswitch:

    2593. 

  2593.             numjl = self.numjl

    2594. 

  2594. 

    2595. 

  2595.         pp = self.P

    2596. 

  2596.         ff = self.F

    2597. 

  2597. 

    2598. 

  2598. 

    2599. 

  2599.         sposli=[]

    2600. 

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

    2601. 

  2601.         # Check if names are valid

    2602. 

  2602.         for name in sntup:

    2603. 

  2603.             if name not in exoli:

    2604. 

  2604.                 return 'Error: '+name+' is not a valid exoshock name for this model!'

    2605. 

  2605.         for name in sntup:

    2606. 

  2606.             sposli.append(exoli.index(name))

    2607. 

  2607.         # Expose spos choice to self for show_irf

    2608. 

  2608.         self.spos = (sntup,sposli)

    2609. 

  2609. 

    2610. 

  2610.         shock = MAT.zeros((nexo,1))

    2611. 

  2611.         sendo = MAT.zeros((nendo,1))

    2612. 

  2612. 

    2613. 

  2613.         for i1,spos in enumerate(sposli):

    2614. 

  2614.             shock[spos,0] = 1.0*shockmod[i1]

    2615. 

  2615.         shock = MAT.vstack((shock,sendo))

    2616. 

  2616.         '''

    2617. 

  2617.         if shockvec == None:

    2618. 

  2618.             for elem in range(ssigma.shape[0]):

    2619. 

  2619.                 ssigma[elem,elem] = np.sqrt(ssigma[elem,elem])

    2620. 

  2620.             shock = np.dot(ssigma,shock)

    2621. 

  2621.         else:

    2622. 

  2622.             ssigma = self.mkssigma(self.tstates,self.nexo,shockvec)

    2623. 

  2623.             for elem in range(ssigma.shape[0]):

    2624. 

  2624.                 ssigma[elem,elem] = np.sqrt(ssigma[elem,elem])            

    2625. 

  2625.             shock = np.dot(ssigma,shock)

    2626. 

  2626.         '''

    2627. 

  2627.         self.sshock = COP.deepcopy(shock)

    2628. 

  2628. 

    2629. 

  2629.         x_one_m1 = shock

    2630. 

  2630.         y_one_0 = ff*x_one_m1

    2631. 

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

    2632. 

  2632.         x_one_0 = pp*x_one_m1

    2633. 

  2633.         if self.oswitch:

    2634. 

  2634.             nother = self.nother

    2635. 

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

    2636. 

  2636. 

    2637. 

  2637.         x_one_c = COP.deepcopy(x_one_m1)

    2638. 

  2638.         y_one_c = COP.deepcopy(y_one_0)

    2639. 

  2639.         if self.oswitch:

    2640. 

  2640.             o_one_c = COP.deepcopy(o_one_0)

    2641. 

  2641. 

    2642. 

  2642.         x_one = COP.deepcopy(x_one_c)

    2643. 

  2643.         y_one = COP.deepcopy(y_one_c)

    2644. 

  2644.         if self.oswitch:

    2645. 

  2645.             o_one = COP.deepcopy(o_one_c)

    2646. 

  2646. 

    2647. 

  2647.         for i1 in range(tlen):

    2648. 

  2648.             x_one_n = pp*x_one_c

    2649. 

  2649.             y_one_n = ff*x_one_c

    2650. 

  2650.             if self.oswitch:

    2651. 

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

    2652. 

  2652. 

    2653. 

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

    2654. 

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

    2655. 

  2655.             if self.oswitch:

    2656. 

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

    2657. 

  2657. 

    2658. 

  2658.             x_one_c = x_one_n

    2659. 

  2659.             y_one_c = y_one_n

    2660. 

  2660. 

    2661. 

  2661.         # Throw away first observation

    2662. 

  2662.         self.irf_x_one = x_one[:,1:-1]

    2663. 

  2663.         self.irf_y_one = y_one[:,1:-1]

    2664. 

  2664.         self.inirf = [self.irf_x_one,self.irf_y_one]

    2665. 

  2665.         if self.oswitch:

    2666. 

  2666.             self.irf_o_one = o_one[:,2:]

    2667. 

  2667.             self.inirf = self.inirf + [self.irf_o_one,]

    2668. 

  2668. 

    2669. 

  2669. #------------------------currently not called and really not working, removed from init----------------

    2670. 

  2670. class FairTaylor:

    2671. 

  2671. 

    2672. 

  2672.     def __init__(self,algtype='type1',intup=None):

    2673. 

  2673.         self.algtype = algtype

    2674. 

  2674.         self.data = intup[0]

    2675. 

  2675.         self.param = intup[1]

    2676. 

  2676.         self.sstates_list = intup[2]

    2677. 

  2677.         self.vardic = intup[3]

    2678. 

  2678.         self.vardic2 = intup[4]

    2679. 

  2679.         self.modeq = intup[5]

    2680. 

  2680.         self.maxlags = intup[6]

    2681. 

  2681.         self.maxleads = intup[7]

    2682. 

  2682.         self.maxshocks = intup[8]

    2683. 

  2683.         self.maxinfosets = intup[9]

    2684. 

  2684.         self.shock_vars = intup[10]

    2685. 

  2685.         self.iid_vars = intup[11]

    2686. 

  2686.         self.vars = intup[12]

    2687. 

  2687.         self.re_var = intup[13]

    2688. 

  2688.         self.re_var2 = intup[14]

    2689. 

  2689.         self.sstate = intup[15]

    2690. 

  2690. 

    2691. 

  2691.         self.prepdat(self.data)

    2692. 

  2692.         self.prepmodeq(self.modeq)

    2693. 

  2693. 

    2694. 

  2694.     def prepdat(self,data=None):

    2695. 

  2695.         self.cur_pres = {}

    2696. 

  2696.         self.cur_past = {}

    2697. 

  2697.         self.cur_fut = {}

    2698. 

  2698.         for x1 in data.items():

    2699. 

  2699.             self.cur_pres[x1[0]] = x1[1][self.maxlags:len(x1[1])-self.maxleads]

    2700. 

  2700.             self.cur_past[x1[0]] = x1[1][0:-2]

    2701. 

  2701.             self.cur_fut[x1[0]] = x1[1][2:]

    2702. 

  2702. 

    2703. 

  2703.     def prepmodeq(self,modeq=None):

    2704. 

  2704.         re_var2 = self.re_var2

    2705. 

  2705.         self.modeq1 = COP.deepcopy(modeq)

    2706. 

  2706.         tmp_list = []

    2707. 

  2707.         i1 = 0

    2708. 

  2708.         for x1 in modeq:

    2709. 

  2709.             itera = re_var2.finditer(x1)

    2710. 

  2710.             for x2 in itera:

    2711. 

  2711.                 if x2.group('svar'):

    2712. 

  2712.                     tmp_list = tmp_list + [[i1,x2.group('svar'),x2.start(),x2.end()],]

    2713. 

  2713.             i1 = i1 + 1

    2714. 

  2714.         tmp_list.reverse()

    2715. 

  2715.         if tmp_list:

    2716. 

  2716.             for x2 in tmp_list:

    2717. 

  2717.                 i1 = x2[0]

    2718. 

  2718.                 expr = x2[1]

    2719. 

  2719.                 estart = x2[2]

    2720. 

  2720.                 eend = x2[3]

    2721. 

  2721.                 x1 = self.modeq1[i1][0:estart]+expr.capitalize()+\

    2722. 

  2722.                    '_bar'+self.modeq1[i1][eend:]

    2723. 

  2723.                 self.modeq1[i1] = x1

    2724. 

  2724.         self.modeq2 = self.modeq1[0:len(modeq)-self.maxshocks]

    2725. 

  2725. 

    2726. 

  2726. 

    2727. 

  2727.         self.modeq3 = COP.deepcopy(self.modeq2)

    2728. 

  2728.         # Create varlist from vardic to establish fixed order

    2729. 

  2729.         self.varlist = self.vardic.items()

    2730. 

  2730.         # Create ordered vars and varnames

    2731. 

  2731.         varnames = []

    2732. 

  2732.         variables = []

    2733. 

  2733.         sub_list1 = []

    2734. 

  2734.         sub_list2 = []

    2735. 

  2735.         tmp_index = {}

    2736. 

  2736.         for x1 in self.varlist:

    2737. 

  2737.             varnames = varnames + [x1[1],]

    2738. 

  2738.             variables = variables + [x1[0],]

    2739. 

  2739.         self.varnames = varnames

    2740. 

  2740.         self.variables = variables

    2741. 

  2741.         re_var2 = self.re_var2

    2742. 

  2742.         i1 = 0

    2743. 

  2743.         i2 = 0

    2744. 

  2744.         for x1 in self.modeq2:

    2745. 

  2745.             itera = re_var2.finditer(x1)

    2746. 

  2746.             for x2 in itera:

    2747. 

  2747.                 if x2.group('vvar') and not x2.group('exp') and not x2.group('lagt'):

    2748. 

  2748.                     tmp_match = variables.index(x2.group('var'))

    2749. 

  2749.                     if not tmp_index.has_key(x2.group('var')):

    2750. 

  2750.                         tmp_index[x2.group('var')] = (i2,x2.group(0))

    2751. 

  2751.                         i2 = i2 + 1

    2752. 

  2752.                     sub_list1 = sub_list1 + [['invar['+str(tmp_index[x2.group('var')][0])+']',

    2753. 

  2753.                                   x2.start(),x2.end(),x2.group(0),tmp_match,i1],]

    2754. 

  2754. 

    2755. 

  2755.             sub_list1.reverse()

    2756. 

  2756.             if sub_list1:

    2757. 

  2757.                 for x4 in sub_list1:

    2758. 

  2758.                     str_tmp = x4[0]

    2759. 

  2759.                     pstart = x4[1]

    2760. 

  2760.                     pend = x4[2]

    2761. 

  2761.                     self.modeq3[i1] = self.modeq3[i1][0:pstart]+str_tmp+self.modeq3[i1][pend:]

    2762. 

  2762.                 sub_list2 = sub_list2 + [sub_list1,]

    2763. 

  2763.                 sub_list1 = []

    2764. 

  2764.             i1 = i1 + 1

    2765. 

  2765.         self.index = tmp_index

    2766. 

  2766.         # Create inverted index2

    2767. 

  2767.         index2 = {}

    2768. 

  2768.         for x1 in self.index.items():

    2769. 

  2769.             index2[x1[1][0]] = (x1[0],x1[1][1])

    2770. 

  2770.         self.index2 = index2

    2771. 

  2771.         return sub_list2

    2772. 

  2772. 

    2773. 

  2773.     def makecur(self,tpos,curlag,curpres,curfut):

    2774. 

  2774.         lag_dic={}

    2775. 

  2775.         pres_dic={}

    2776. 

  2776.         fut_dic={}

    2777. 

  2777.         all_dic={}

    2778. 

  2778.         for x1 in curlag.items():

    2779. 

  2779.             if float(x1[1][tpos].data) > 0:

    2780. 

  2780.                 lag_dic[self.vardic2[x1[0]]+'_L'] = float(x1[1][tpos].data)

    2781. 

  2781.             else:

    2782. 

  2782.                 lag_dic[self.vardic2[x1[0]]+'_L'] = 0.01

    2783. 

  2783.         for x1 in curpres.items():

    2784. 

  2784.             if float(x1[1][tpos].data) > 0:

    2785. 

  2785.                 pres_dic[self.vardic2[x1[0]]+'_N'] = float(x1[1][tpos].data)

    2786. 

  2786.             else:

    2787. 

  2787.                 pres_dic[self.vardic2[x1[0]]+'_N'] = 0.01

    2788. 

  2788.         for x1 in curfut.items():

    2789. 

  2789.             if float(x1[1][tpos].data) > 0:

    2790. 

  2790.                 fut_dic['INFO_N_'+self.vardic2[x1[0]]+'_F'] = float(x1[1][tpos].data)

    2791. 

  2791.             else:

    2792. 

  2792.                 fut_dic['INFO_N_'+self.vardic2[x1[0]]+'_F'] = 0.01

    2793. 

  2793.         all_dic.update(lag_dic)

    2794. 

  2794.         all_dic.update(pres_dic)

    2795. 

  2795.         all_dic.update(fut_dic)

    2796. 

  2796.         self.alldic =all_dic

    2797. 

  2797.         return all_dic

    2798. 

  2798. 

    2799. 

  2799.     def solveone(self,tpos,curlag,curpres,curfut):

    2800. 

  2800. 

    2801. 

  2801.         #Define the function to be handed over

    2802. 

  2802.         #to fsolve

    2803. 

  2803. 

    2804. 

  2804.         def func(invar):

    2805. 

  2805.             locals().update(self.param)

    2806. 

  2806.             locals().update(self.sstate[0])

    2807. 

  2807.             locals().update(self.makecur(tpos,self.cur_past,self.cur_pres,self.cur_fut))

    2808. 

  2808.             fdot1 = np.zeros((len(self.modeq3)),'d')

    2809. 

  2809.             i1=0

    2810. 

  2810.             for x in self.modeq3:

    2811. 

  2811.                 fdot1[i1] = eval(x)

    2812. 

  2812.                 i1 = i1 + 1

    2813. 

  2813.             return fdot1

    2814. 

  2814. 

    2815. 

  2815. 

    2816. 

  2816.         #def func(invar):

    2817. 

  2817.             #locals().update(self.param)

    2818. 

  2818.             #locals().update(self.sstate[0])

    2819. 

  2819.             #locals().update(self.makecur(tpos,self.cur_past,self.cur_pres,self.cur_fut))

    2820. 

  2820.             #fdot1 = np.zeros((len(self.modeq3)),'d')

    2821. 

  2821.             #i1=0

    2822. 

  2822.             #for x in self.modeq3:

    2823. 

  2823.                 #fdot1[i1] = abs(eval(x))

    2824. 

  2824.                 #i1 = i1 + 1

    2825. 

  2825.             #return sum(fdot1)

    2826. 

  2826. 

    2827. 

  2827.         #def fprime(invar):

    2828. 

  2828.             #return np.array(([1,-1]))

    2829. 

  2829. 

    2830. 

  2830. 

    2831. 

  2831.         # Define the initial values and

    2832. 

  2832.         # start the non-linear solver

    2833. 

  2833.         init_val = len(self.index)*[0,]

    2834. 

  2834.         for x1 in self.index.items():

    2835. 

  2835.             varia = self.vardic[x1[0]]['var']

    2836. 

  2836.             init_val[x1[1][0]] = float(curpres[varia][tpos].data)

    2837. 

  2837. 

    2838. 

  2838.         # Determine bounds

    2839. 

  2839.         #in_bounds = len(self.index)*[(0.01,None),]

    2840. 

  2840. 

    2841. 

  2841.         (output,infodict,ier,mesg) = O.fsolve(func,init_val,full_output=1)

    2842. 

  2842. 

    2843. 

  2843.         #(output,f_out,d_out) = O.fmin_l_bfgs_b(func,init_val,

    2844. 

  2844.                                                             #fprime,approx_grad=True,bounds=in_bounds)

    2845. 

  2845. 

    2846. 

  2846. 

    2847. 

  2847.         # Attach the outputs of the solver as attributes

    2848. 

  2848.         fsout={}

    2849. 

  2849.         self.output = output

    2850. 

  2850.         i1 = 0

    2851. 

  2851.         for x1 in self.index2.items():

    2852. 

  2852.             fsout[x1[1][1]] = (output[x1[0]],x1[1][0])

    2853. 

  2853.             i1 = i1 + 1

    2854. 

  2854.         self.fsout = fsout

    2855. 

  2855.         self.output = output

    2856. 

  2856.         self.infodict = infodict

    2857. 

  2857.         self.ier = ier

    2858. 

  2858.         self.mesg = mesg

    2859. 

  2859. 

    2860. 

  2860.         return fsout

    2861. 

  2861. 

    2862. 

  2862.     def solveall(self):

    2863. 

  2863.         # Calculate initial convergence criterion

    2864. 

  2864.         criterion = 10

    2865. 

  2865.         loop_c = 0

    2866. 

  2866.         while criterion > 0.01:

    2867. 

  2867.             cur_past_t = {}

    2868. 

  2868.             cur_pres_t = {}

    2869. 

  2869.             cur_fut_t = {}

    2870. 

  2870.             for x1 in self.vardic2.items():

    2871. 

  2871.                 cur_past_t[x1[0]] = np.zeros(len(self.cur_pres.items()[0][1]))

    2872. 

  2872.                 cur_past_t[x1[0]][0] = float(self.cur_past[x1[0]][0].data)

    2873. 

  2873.                 cur_pres_t[x1[0]] = np.zeros(len(self.cur_pres.items()[0][1]))

    2874. 

  2874.                 cur_fut_t[x1[0]] = np.zeros(len(self.cur_pres.items()[0][1]))

    2875. 

  2875.                 cur_fut_t[x1[0]][-1] = float(self.cur_fut[x1[0]][-1].data)

    2876. 

  2876.             for i1 in range(0,len(self.cur_pres.items()[0][1])):

    2877. 

  2877.                 out_tmp = self.solveone(i1,self.cur_past,self.cur_pres,self.cur_fut)

    2878. 

  2878.                 for x2 in out_tmp.items():

    2879. 

  2879.                     # Standard

    2880. 

  2880.                     if self.algtype == 'type1':

    2881. 

  2881.                         cur_pres_t[self.vardic[x2[1][1]['var']]][i1] = x2[1][0]

    2882. 

  2882.                     # Fast-Gauss Seidel

    2883. 

  2883.                     elif self.algtype == 'type2':

    2884. 

  2884.                         cur_pres_t[self.vardic[x2[1][1]['var']]][i1] = x2[1][0]

    2885. 

  2885.                         try:

    2886. 

  2886.                             self.cur_past[self.vardic[x2[1][1]['var']]][i1+1] = x2[1][0]

    2887. 

  2887.                         except:

    2888. 

  2888.                             pass

    2889. 

  2889. 

    2890. 

  2890.             # Recalculate convergence criterion

    2891. 

  2891.             cr_list = [0,]*len(self.vardic)

    2892. 

  2892.             va_list = [0,]*len(self.vardic)

    2893. 

  2893.             for x1,x2 in zip(self.vardic.items(),range(0,len(self.vardic),1)):

    2894. 

  2894.                 cr_list[x2] = cr_list[x2] + sum(abs(cur_fut_t[x1[1]][0:-1]-cur_pres_t[x1[1]][1:]))

    2895. 

  2895.                 va_list[x2] = x1[1]

    2896. 

  2896.             criterion = sum(cr_list)

    2897. 

  2897.             loop_c = loop_c + 1

    2898. 

  2898.             print ('Loop: ',loop_c)

    2899. 

  2899.             print ('Var:  ',va_list)

    2900. 

  2900.             print ('Crit: ',cr_list)

    2901. 

  2901. 

    2902. 

  2902. 

    2903. 

  2903.             for x1 in self.vardic2.items():

    2904. 

  2904.                 cur_past_t[x1[0]][1:] = cur_pres_t[x1[0]][:-1]

    2905. 

  2905.                 cur_fut_t[x1[0]][0:-1] = cur_pres_t[x1[0]][1:]

    2906. 

  2906. 

    2907. 

  2907.             for x1,x2,x3 in zip(cur_past_t.items(),cur_pres_t.items(),cur_fut_t.items()):

    2908. 

  2908.                 self.cur_past[x1[0]] = TSS.time_series\

    2909. 

  2909.                     (x1[1],start_date=TSS.Date(freq=self.cur_past[x1[0]].start_date.freqstr,

    2910. 

  2910.                                    year=self.cur_past[x1[0]].start_date.year,

    2911. 

  2911.                                    quarter=self.cur_past[x1[0]].start_date.quarter),

    2912. 

  2912.                                    freq=self.cur_past[x1[0]].start_date.freqstr)

    2913. 

  2913. 

    2914. 

  2914.                 self.cur_pres[x1[0]] = TSS.time_series\

    2915. 

  2915.                     (x2[1],start_date=TSS.Date(freq=self.cur_pres[x1[0]].start_date.freqstr,

    2916. 

  2916.                                    year=self.cur_pres[x1[0]].start_date.year,

    2917. 

  2917.                                    quarter=self.cur_pres[x1[0]].start_date.quarter),

    2918. 

  2918.                                    freq=self.cur_pres[x1[0]].start_date.freqstr)

    2919. 

  2919. 

    2920. 

  2920.                 self.cur_fut[x1[0]] = TSS.time_series\

    2921. 

  2921.                     (x3[1],start_date=TSS.Date(freq=self.cur_fut[x1[0]].start_date.freqstr,

    2922. 

  2922.                                    year=self.cur_fut[x1[0]].start_date.year,

    2923. 

  2923.                                    quarter=self.cur_fut[x1[0]].start_date.quarter),

    2924. 

  2924.                                    freq=self.cur_fut[x1[0]].start_date.freqstr)

    2925. 

  2925.