In [1]:

    %matplotlib inline

    import numpy as np
    import sympy as sy
    import pandas as pd
    import statsmodels.api as sm
    import statsmodels.formula.api as smf
    import matplotlib.pyplot as plt

    from numpy.linalg import eig, inv
    from pandas_datareader.data import DataReader

    sy.init_printing(use_latex='mathjax')
    np.set_printoptions(precision=4, suppress=True)

    C:\Users\Marcin\Anaconda3\lib\site-packages\statsmodels\compat\pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.
      from pandas.core import datetools

In [2]:

    from matplotlib import rcParams

    # Restore old behavior of rounding default axis ranges
    rcParams['axes.autolimit_mode'] = 'round_numbers'
    rcParams['axes.xmargin'] = 0
    rcParams['axes.ymargin'] = 0

    # Adjust tick placement
    rcParams['xtick.direction'] = 'in'
    rcParams['ytick.direction'] = 'in'
    rcParams['xtick.top'] = True
    rcParams['ytick.right'] = True

    # Disable legend frame
    rcParams['legend.frameon'] = False

In [3]:

    rec = DataReader('USREC', 'fred', start='1947-01', end='2020-01')
    rec = rec.dropna()

    pop = DataReader(['B230RC0Q173SBEA','CNP16OV'], 'fred', start='1947-01', end='2020-01')
    pop = pop.dropna()
    pop = pop.resample('QS').mean()
    pop = pop.dropna()

    pop.tail()

Out[3]:

| B230RC0Q173SBEA | CNP16OV
--- | --- | ---
DATE |  |
2017-01-01 | 325108.0 | 254082.0
2017-04-01 | 325640.0 | 254588.0
2017-07-01 | 326276.0 | 255151.0
2017-10-01 | 326907.0 | 255766.0
2018-01-01 | 327423.0 | 256780.0

In [4]:

    fred = DataReader(['PRS85006053','GDP','PCND','PCESV','PCDG',
                       'FPI','HOANBS','COMPNFB','GDPDEF','TB3MS'],
                      'fred', start='1947-01', end='2020-01')
    fred = fred.dropna()
    fred = fred.resample('QS').mean()
    fred = fred.dropna()
    fred.head()

Out[4]:

| PRS85006053 | GDP | PCND | PCESV | PCDG | FPI | HOANBS | COMPNFB | GDPDEF | TB3MS
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
DATE |  |  |  |  |  |  |  |  |  |
1947-01-01 | 1.696 | 243.080 | 74.872 | 60.695 | 20.721 | 35.361 | 48.186 | 3.971 | 12.566 | 0.38
1947-04-01 | 1.759 | 246.267 | 76.898 | 61.908 | 21.352 | 35.744 | 48.190 | 4.062 | 12.745 | 0.38
1947-07-01 | 1.787 | 250.115 | 78.648 | 63.249 | 21.770 | 37.826 | 48.320 | 4.181 | 12.957 | 0.66
1947-10-01 | 1.873 | 260.309 | 79.966 | 64.344 | 23.487 | 41.793 | 48.792 | 4.280 | 13.276 | 0.85
1948-01-01 | 1.927 | 266.173 | 81.546 | 65.406 | 23.548 | 43.574 | 49.103 | 4.382 | 13.379 | 0.97

In [5]:

    fred.tail()

Out[5]:

| PRS85006053 | GDP | PCND | PCESV | PCDG | FPI | HOANBS | COMPNFB | GDPDEF | TB3MS
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
DATE |  |  |  |  |  |  |  |  |  |
2017-01-01 | 136.035 | 19057.705 | 2787.636 | 8960.721 | 1443.205 | 3128.877 | 113.252 | 117.102 | 112.746 | 0.51
2017-04-01 | 137.560 | 19250.009 | 2790.618 | 9059.818 | 1456.586 | 3173.297 | 113.864 | 117.254 | 113.029 | 0.80
2017-07-01 | 139.618 | 19500.602 | 2823.815 | 9127.685 | 1477.617 | 3207.271 | 114.240 | 118.300 | 113.614 | 1.07
2017-10-01 | 141.730 | 19754.102 | 2883.975 | 9252.544 | 1517.802 | 3279.209 | 115.178 | 119.018 | 114.275 | 1.07
2018-01-01 | 143.311 | 19965.326 | 2912.268 | 9370.815 | 1499.175 | 3340.484 | 115.770 | 120.007 | 114.837 | 1.41

In [6]:

    pop_hp_cycle, pop_hp_trend = sm.tsa.filters.hpfilter(np.log(pop), lamb=10000)
    pop_smooth = np.exp(pop_hp_trend)

    fig, ax = plt.subplots()

    pop['CNP16OV'].to_period('D').plot(ax=ax, lw=2)
    pop_smooth['CNP16OV'].to_period('D').plot(ax=ax, lw=2)
    plt.legend(['raw', 'filtered'])
    plt.show()

    fig, ax = plt.subplots()

    pop['CNP16OV'].pct_change().to_period('D').plot(ax=ax, lw=2)
    pop_smooth['CNP16OV'].pct_change().to_period('D').plot(ax=ax, lw=2)
    plt.legend(['raw', 'filtered'], loc='upper left')
    plt.show()

In [7]:

    dta = fred[['GDP']]
    dta.is_copy = False

    dta['Output'] = np.log(fred['GDP']*10**9/fred['GDPDEF']*100
                           /(pop_smooth['CNP16OV']*10**3))
    dta['Consumption'] = np.log((fred['PCND']+fred['PCESV'])*10**9/fred['GDPDEF']*100
                                /(pop_smooth['CNP16OV']*10**3))
    dta['Investment'] = np.log((fred['PCDG']+fred['FPI'])*10**9/fred['GDPDEF']*100
                               /(pop_smooth['CNP16OV']*10**3))
    dta['Capital'] = 0*dta['Output']
    dta['Hours'] = np.log(fred['HOANBS']*100*fred['GDP']/np.mean(fred['GDP']['2010-01':'2010-10'])
                          /fred['PRS85006053']/pop_smooth['CNP16OV'])
    dta['Wages'] = np.log(fred['COMPNFB']/fred['GDPDEF']*100)
    dta['Interest Rate'] = (1+fred['TB3MS']/100)**(1/4)/(1+fred['GDPDEF'].pct_change())
    dta['TFP'] = 0*dta['Output']
    dta['Productivity'] = np.log(fred['PRS85006053']/fred['GDPDEF']*100)-dta['Hours']
    dta['Price Level'] = np.log(fred['GDPDEF'])

    dta = dta.drop('GDP', 1)
    dta = dta.dropna()

    dta.head()

Out[7]:

| Output | Consumption | Investment | Capital | Hours | Wages | Interest Rate | TFP | Productivity | Price Level
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
DATE |  |  |  |  |  |  |  |  |  |
1948-01-01 | 9.874765 | 9.280724 | 8.497130 | 0.0 | -7.722532 | 3.488989 | 0.994699 | 0.0 | 10.389980 | 2.593686
1948-04-01 | 9.888965 | 9.292191 | 8.506602 | 0.0 | -7.713483 | 3.493807 | 0.993726 | 0.0 | 10.385552 | 2.602467
1948-07-01 | 9.892542 | 9.282950 | 8.517118 | 0.0 | -7.711741 | 3.495513 | 0.984260 | 0.0 | 10.392732 | 2.620821
1948-10-01 | 9.891661 | 9.285558 | 8.507375 | 0.0 | -7.727315 | 3.502124 | 0.999734 | 0.0 | 10.413690 | 2.623871
1949-01-01 | 9.875913 | 9.285166 | 8.465233 | 0.0 | -7.751860 | 3.513677 | 1.008177 | 0.0 | 10.430040 | 2.618636

In [8]:

    dta.tail()

Out[8]:

| Output | Consumption | Investment | Capital | Hours | Wages | Interest Rate | TFP | Productivity | Price Level
--- | --- | --- | --- | --- | --- | --- | --- | --- | --- | ---
DATE |  |  |  |  |  |  |  |  |  |
2017-01-01 | 11.104037 | 10.620279 | 9.676534 | 0.0 | -7.783061 | 4.643078 | 0.996335 | 0.0 | 12.576006 | 4.725138
2017-04-01 | 11.109106 | 10.623959 | 9.684126 | 0.0 | -7.781244 | 4.641868 | 0.999485 | 0.0 | 12.582830 | 4.727644
2017-07-01 | 11.114414 | 10.624826 | 9.688311 | 0.0 | -7.782327 | 4.645587 | 0.997502 | 0.0 | 12.593600 | 4.732807
2017-10-01 | 11.119066 | 10.631924 | 9.703698 | 0.0 | -7.778710 | 4.645837 | 0.996865 | 0.0 | 12.599196 | 4.738608
2018-01-01 | 11.122333 | 10.636559 | 9.705180 | 0.0 | -7.776503 | 4.649207 | 0.998595 | 0.0 | 12.603177 | 4.743514

In [9]:

    # Estimate Capital series using PIM
    temp = fred[['GDP']]
    temp.is_copy = False

    temp['Inv'] = (fred['PCDG']+fred['FPI'])/fred['GDPDEF']*100/4
    temp['LnInv'] = np.log(temp['Inv'])
    temp['t'] = np.arange(len(temp['Inv']))

    trend = smf.ols(formula='LnInv ~ t', data=temp).fit()
    intercept, slope = trend.params

    delta = 0.025
    K = np.zeros(len(temp['LnInv']))
    K_init = np.exp(intercept-slope)/(slope+delta)
    K[0] = (1-delta)*K_init+temp['Inv'][0]
    for i in range(1,len(temp['Inv'])):
        K[i] = (1-delta)*K[i-1]+temp['Inv'][i]
    temp['Cap'] = K

    temp['Cap'].to_period('D').plot(lw=2)
    plt.show()

    (temp['Cap']/(fred['GDP']/fred['GDPDEF']*100)).to_period('D').plot(lw=2)
    plt.show()

In [10]:

    dta['Capital'] = np.log(temp['Cap']*10**9/(pop_smooth['CNP16OV']*10**3))

    α = 1/3
    dta['TFP'] = dta['Output']-α*dta['Capital']-(1-α)*dta['Hours']
    dta['TFP'].to_period('D').plot(lw=2)
    plt.show()

In [11]:

    RGDP_pc = (fred['GDP']*10**9/fred['GDPDEF']*100/(pop_smooth['B230RC0Q173SBEA']*10**3)).dropna()

    fig, ax = plt.subplots()

    RGDP_pc.to_period('D').plot(ax=ax, lw=2, style='k-')

    ax.set_ylim(10000, 60000)
    ylim = ax.get_ylim()

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP (2009 dollars)')
    plt.show()

In [12]:

    fig, ax = plt.subplots()

    np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')

    ticks = [10000, 20000, 30000, 40000, 50000, 60000]
    ax.set_yticks(np.log(ticks))
    ax.set_yticklabels(ticks)

    ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))

    ylim = ax.get_ylim()
    ax.set_ylim(ylim)

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, log scale')
    plt.show()

In [13]:

    fig, ax = plt.subplots()

    (100*RGDP_pc.pct_change(4)).to_period('D').plot(ax=ax, lw=2, style='k-')

    avg = np.mean(100*RGDP_pc.to_period('D').pct_change(4))

    ax.set_ylim(-6, 12)
    ylim = ax.get_ylim()

    ax.hlines(avg, dta.index[0], dta.index[-1], color='r', linewidth=2)
    ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, year-over-year change (%)')
    plt.show()

    print('Average growth rate (%) =', avg)

    Average growth rate (%) = 1.987401422162432

In [14]:

    gdp_lin_cycle, gdp_lin_trend = sm.tsa.filters.hpfilter(np.log(RGDP_pc), lamb=1e9)
    gdp_hp_cycle, gdp_hp_trend = sm.tsa.filters.hpfilter(np.log(RGDP_pc))
    gdp_cf_cycle, gdp_cf_trend = sm.tsa.filters.cffilter(np.log(RGDP_pc))

In [15]:

    fig, ax = plt.subplots()

    np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')
    gdp_lin_trend.plot(ax=ax, lw=2, style='r-')

    ticks = [10000, 20000, 30000, 40000, 50000, 60000]
    ax.set_yticks(np.log(ticks))
    ax.set_yticklabels(ticks)

    ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))

    ylim = ax.get_ylim()
    ax.set_ylim(ylim)

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, exponential trend')
    plt.show()

In [16]:

    fig, ax = plt.subplots()

    (100*gdp_lin_cycle).to_period('D').plot(ax=ax, lw=2, style='k-')
    ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)

    ylim = ax.get_ylim()
    ax.set_ylim(ylim)

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, exponential trend residual (%)')
    plt.show()

In [17]:

    fig, ax = plt.subplots()

    np.log(RGDP_pc).to_period('D').plot(ax=ax, lw=2, style='k-')
    gdp_hp_trend.plot(ax=ax, lw=2, style='r-')

    ticks = [10000, 20000, 30000, 40000, 50000, 60000]
    ax.set_yticks(np.log(ticks))
    ax.set_yticklabels(ticks)

    ax.set_ylim(np.log(ticks[0]), np.log(ticks[-1]))

    ylim = ax.get_ylim()
    ax.set_ylim(ylim)

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, Hodrick-Prescott trend')
    plt.show()

In [18]:

    fig, ax = plt.subplots()

    (100*gdp_hp_cycle.to_period('D')).plot(ax=ax, lw=2, style='k-')
    ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)

    ax.set_ylim(-6, 6)
    ylim = ax.get_ylim()

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, Hodrick-Prescott residual (%)')
    plt.show()

In [19]:

    fig, ax = plt.subplots()

    (100*gdp_cf_cycle).to_period('D').plot(ax=ax, lw=2, style='k-')
    ax.hlines(0, dta.index[0], dta.index[-1], linewidth=0.5)

    ylim = ax.get_ylim()

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.title('US real per capita GDP, Christiano-Fitzgerald residual (%)')
    plt.show()

In [20]:

    np.corrcoef(gdp_hp_cycle, gdp_cf_cycle)

Out[20]:

    array([[1.    , 0.8881],
           [0.8881, 1.    ]])

In [21]:

    hp_cycles, hp_trend = sm.tsa.filters.hpfilter((100*dta).dropna())
    cf_cycles, cf_trend = sm.tsa.filters.cffilter((100*dta).dropna())

In [22]:

    print('Standard Deviations')
    print(hp_cycles.std())

    print('')
    print('Autocorrelations')
    a = list(dta.columns.values)
    for i in range(len(a)):
        print(dta.columns.values[i], '  \t\t', hp_cycles[dta.columns.values[i]].autocorr())

    print('')
    print('Correlations')
    print(hp_cycles.corr(method='pearson'))

    Standard Deviations
    Output           1.602339
    Consumption      0.858648
    Investment       4.539281
    Capital          0.571145
    Hours            1.599009
    Wages            0.834119
    Interest Rate    0.392727
    TFP              0.996028
    Productivity     1.302317
    Price Level      0.888267
    dtype: float64

    Autocorrelations
    Output   		 0.8497228502803469
    Consumption   		 0.8319719570655278
    Investment   		 0.8686667132636827
    Capital   		 0.966698644398374
    Hours   		 0.9009160329580754
    Wages   		 0.6546991130083243
    Interest Rate   		 0.397955895107275
    TFP   		 0.7052710432683185
    Productivity   		 0.6470773019467758
    Price Level   		 0.9082135152943024

    Correlations
                     Output  Consumption  Investment   Capital     Hours  \
    Output         1.000000     0.752088    0.791239  0.364096  0.810409
    Consumption    0.752088     1.000000    0.641177  0.397785  0.685424
    Investment     0.791239     0.641177    1.000000  0.242497  0.634787
    Capital        0.364096     0.397785    0.242497  1.000000  0.597770
    Hours          0.810409     0.685424    0.634787  0.597770  1.000000
    Wages          0.095967     0.212493    0.053883  0.099806 -0.091599
    Interest Rate -0.013552     0.035548   -0.068175  0.037618  0.013819
    TFP            0.671789     0.400292    0.547151 -0.245177  0.119213
    Productivity   0.506890     0.302061    0.529303 -0.232276 -0.060970
    Price Level   -0.146255    -0.250324   -0.393407  0.175168  0.084581

                      Wages  Interest Rate       TFP  Productivity  Price Level
    Output         0.095967      -0.013552  0.671789      0.506890    -0.146255
    Consumption    0.212493       0.035548  0.400292      0.302061    -0.250324
    Investment     0.053883      -0.068175  0.547151      0.529303    -0.393407
    Capital        0.099806       0.037618 -0.245177     -0.232276     0.175168
    Hours         -0.091599       0.013819  0.119213     -0.060970     0.084581
    Wages          1.000000      -0.002761  0.233342      0.253733    -0.248189
    Interest Rate -0.002761       1.000000 -0.043781      0.007753    -0.016279
    TFP            0.233342      -0.043781  1.000000      0.925099    -0.359290
    Productivity   0.253733       0.007753  0.925099      1.000000    -0.393792
    Price Level   -0.248189      -0.016279 -0.359290     -0.393792     1.000000

In [23]:

    cf_cycles['Investment / 4'] = cf_cycles['Investment'] / 4

    fig, ((ax1, ax2, ax3), (ax4, ax5, ax6),
          (ax7, ax8, ax9)) = plt.subplots(3, 3, figsize=(20,15), sharex=False, sharey=False)

    cf_cycles[['Output','Consumption']].to_period('D').plot(ax=ax1, style=['k','r'])

    cf_cycles[['Output','Investment / 4']].to_period('D').plot(ax=ax2, style=['k','r'])

    cf_cycles[['Output','Capital']].to_period('D').plot(ax=ax3, style=['k','r'])

    cf_cycles[['Output','Hours']].to_period('D').plot(ax=ax4, style=['k','r'])

    cf_cycles[['Output','Wages']].to_period('D').plot(ax=ax5, style=['k','r'])

    cf_cycles[['Output','Interest Rate']].to_period('D').plot(ax=ax6, style=['k','r'])

    cf_cycles[['Output','TFP']].to_period('D').plot(ax=ax7, style=['k','r'])

    cf_cycles[['Output','Productivity']].to_period('D').plot(ax=ax8, style=['k','r'])

    cf_cycles[['Output','Price Level']].to_period('D').plot(ax=ax9, style=['k','r'])

    plt.savefig('US_CF.pdf', transparent=True, bbox_inches='tight', pad_inches=0.05)
    plt.show()

In [24]:

    fig, ax = plt.subplots()

    cf_cycles[['Output','Consumption','Investment / 4','Hours']].to_period('D').plot(ax=ax)

    ylim = ax.get_ylim()

    ax.fill_between(rec.index, ylim[0], ylim[1], rec['USREC'], facecolor='lightgrey', edgecolor='lightgrey')

    plt.legend(ncol=2, frameon=True)
    plt.show()

In [25]:

    dta['t'] = np.arange(len(dta['TFP']))
    trend_TFP = smf.ols(formula='TFP ~ t', data=dta).fit()
    dta['TFP_resid'] = trend_TFP.resid
    intercept_TFP, slope_TFP = trend_TFP.params
    print(trend_TFP.summary())

                                OLS Regression Results
    ==============================================================================
    Dep. Variable:                    TFP   R-squared:                       0.971
    Model:                            OLS   Adj. R-squared:                  0.971
    Method:                 Least Squares   F-statistic:                     9320.
    Date:                Tue, 08 May 2018   Prob (F-statistic):          2.11e-216
    Time:                        00:36:45   Log-Likelihood:                 505.54
    No. Observations:                 281   AIC:                            -1007.
    Df Residuals:                     279   BIC:                            -999.8
    Df Model:                           1
    Covariance Type:            nonrobust
    ==============================================================================
                     coef    std err          t      P>|t|      [0.025      0.975]
    ------------------------------------------------------------------------------
    Intercept     11.6098      0.005   2428.394      0.000      11.600      11.619
    t              0.0029   2.95e-05     96.542      0.000       0.003       0.003
    ==============================================================================
    Omnibus:                       17.396   Durbin-Watson:                   0.039
    Prob(Omnibus):                  0.000   Jarque-Bera (JB):               29.102
    Skew:                          -0.379   Prob(JB):                     4.79e-07
    Kurtosis:                       4.382   Cond. No.                         323.
    ==============================================================================

    Warnings:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

In [26]:

    dta['TFP'].to_period('D').plot(lw=2)
    (intercept_TFP + slope_TFP*dta['t']).plot(lw=2)
    plt.show()

In [27]:

    dta['TFP_resid'].to_period('D').plot(lw=2)
    plt.show()

In [28]:

    resid_TFP = smf.ols(formula='TFP_resid ~ TFP_resid.shift() -1', data=dta).fit()
    print(resid_TFP.summary())

    print('')
    print('TFP residual autocorrelation    =', resid_TFP.params[0])
    print('TFP residual standard deviation =', resid_TFP.resid.std())

                                OLS Regression Results
    ==============================================================================
    Dep. Variable:              TFP_resid   R-squared:                       0.960
    Model:                            OLS   Adj. R-squared:                  0.960
    Method:                 Least Squares   F-statistic:                     6766.
    Date:                Tue, 08 May 2018   Prob (F-statistic):          1.18e-197
    Time:                        00:36:45   Log-Likelihood:                 960.95
    No. Observations:                 280   AIC:                            -1920.
    Df Residuals:                     279   BIC:                            -1916.
    Df Model:                           1
    Covariance Type:            nonrobust
    =====================================================================================
                            coef    std err          t      P>|t|      [0.025      0.975]
    -------------------------------------------------------------------------------------
    TFP_resid.shift()     0.9622      0.012     82.254      0.000       0.939       0.985
    ==============================================================================
    Omnibus:                        6.450   Durbin-Watson:                   1.925
    Prob(Omnibus):                  0.040   Jarque-Bera (JB):                9.515
    Skew:                           0.091   Prob(JB):                      0.00859
    Kurtosis:                       3.885   Cond. No.                         1.00
    ==============================================================================

    Warnings:
    [1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

    TFP residual autocorrelation    = 0.9621590594231401
    TFP residual standard deviation = 0.007828026718978089

In [29]:

    class Dynare(object):

        def __init__(self, var, varexo, param_values, model, initval):
            self.var = var
            self.varexo = varexo
            self.param_values = param_values
            self.model = model
            self.initval = initval

            self.TranslateInputs()

            return None

        def TranslateInputs(self):
            """
            Convert input strings into sympy objects
            and translate timing
            """

            # Endogenous variables
            self.var_symbols = sy.symbols(self.var)

            self.n = len(self.var_symbols)

            self.current = ()
            self.future  = ()
            self.past    = ()
            self.steadys = ()

            for i in range(self.n):
                self.current += (sy.symbols(str(self.var_symbols[i])+str('_{t}'))),
                self.future  += (sy.symbols(str(self.var_symbols[i])+str('_{t+1}'))),
                self.past    += (sy.symbols(str(self.var_symbols[i])+str('_{t-1}'))),
                self.steadys += (sy.symbols(str(self.var_symbols[i])+str('_{ss}'))),

            # Exogenous variables
            self.varexo_symbols = sy.symbols(self.varexo),

            self.p = len(self.varexo_symbols)

            self.shocks = ()

            for i in range(self.p):
                self.shocks += (sy.symbols(str(self.varexo_symbols[i])+str('_{t}'))),

            # Model equations
            self.symbol_dict = {'betta':sy.symbols('beta'), 'gama':sy.symbols('gamma')}
            self.timing_dict = {}

            for i in range(self.n):
                self.timing_dict[sy.sympify(str(self.var_symbols[i])+'(+1)')] = self.future[i]
                self.timing_dict[sy.sympify(str(self.var_symbols[i])+'(-1)')] = self.past[i]
                self.timing_dict[sy.sympify(str(self.var_symbols[i]))]        = self.current[i]

            for i in range(self.p):
                self.timing_dict[sy.sympify(str(self.varexo_symbols[i]))] = self.shocks[i]

            self.model_symbols = sy.sympify(self.model)
            self.model_symbols = self.model_symbols.subs(self.timing_dict)
            self.model_symbols = self.model_symbols.subs(self.symbol_dict)

            try:
                temp = len(self.model_symbols)
            except:
                self.model_symbols = self.model_symbols,

            self.system = sy.Matrix(self.model_symbols)

            return None

        def SteadySystem(self):
            """
            Removes lead/lag structure from the model to prepare for steady state calculation
            """

            self.steady_vars = {}
            for i in range(self.n):
                self.steady_vars[self.current[i]] = self.steadys[i]
                self.steady_vars[self.future[i]]  = self.steadys[i]
                self.steady_vars[self.past[i]]    = self.steadys[i]
            for i in range(self.p):
                self.steady_vars[self.shocks[i]]  = 0

            self.steady_system = self.system.subs(self.steady_vars)

            return self.steady_system

        def SteadyValues(self):
            """
            Numerically solves for the steady state of the system given initval
            """

            self.SteadySystem()

            try:
                ss = sy.nsolve(self.steady_system.subs(self.param_values), self.steadys, self.initval)
            except:
                raise RuntimeError('Adjust initial values')
            ss = ss.T.tolist()

            self.steady_values = {}
            for i in range(self.n):
                self.steady_values[self.steadys[i]] = np.float(ss[0][i])

            return self.steady_values

        def steady(self):
            """
            Prints out steady state values for the user
            """

            self.SteadyValues()

            print('\n' + 'STEADY-STATE RESULTS' + '\n')

            for i in range(self.n):
                print(str(self.var_symbols[i]), '\t%.4f' % self.steady_values[self.steadys[i]])

            return None

        def resid(self):

            self.SteadyValues()

            temp = self.steady_system.subs(self.param_values).subs(self.steady_values)

            print('\n' + 'Residuals of the static equations' + '\n')

            for i in range(self.n):
                print('Equation number', i, ': %.4f' % temp[i])

            return None

        def TimeIteration(self):
            """
            Solves the first-order approximation of the model using time iteration (thanks to Pontus Rendahl)
            """

            self.SteadyValues()

            self.A_symb = self.system.jacobian(self.past)
            self.B_symb = self.system.jacobian(self.current)
            self.C_symb = self.system.jacobian(self.future)
            self.D_symb = self.system.jacobian(self.shocks)

            self.A = np.array(self.A_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
            self.B = np.array(self.B_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
            self.C = np.array(self.C_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)
            self.D = np.array(self.D_symb.subs(self.steady_vars).subs(self.steady_values).subs(self.param_values)).astype(float)

            self.metric = 1
            self.F = np.zeros((self.n, self.n))
            self.S = np.zeros((self.n, self.n))

            # Add maxit to while loop?

            while self.metric > 1e-13:
                self.F = inv(self.B + self.C @ self.F) @ (-self.A)
                self.S = inv(self.B + self.A @ self.S) @ (-self.C)

                self.metric = np.max(np.max(np.abs(self.A + self.B @ self.F + self.C @ self.F @ self.F)))

            self.Q = -inv(self.B + self.C @ self.F) @ self.D

            # Need formal BK check?

            if sum(eig(self.F)[0] > 1) != 0:
                raise RuntimeError('Blanchard Kahn conditions are not satisfied: no stable equilibrium')
            if sum(eig(self.S)[0] > 1) != 0:
                raise RuntimeError('Blanchard Kahn conditions are not satisfied: indeterminacy')

            return None

        def SimulatedMoments(self, hp_filter=None, shocks_stderr=0.01, periods=10000):

            self.TimeIteration()

            x = np.zeros((self.n, periods))
            ɛ = np.zeros((self.p, periods))

            for i in range(self.p):
                ɛ[i, :] = shocks_stderr * np.random.randn(periods)

            for t in range(1, periods):
                x[:, t] = self.F @ x[:, t-1] + self.Q @ ɛ[:, t]

            print('SIMULATED MOMENTS')
            print('')
            print('VARIABLE \t STD. DEV.')

            if hp_filter == None:
                for i in range(self.n):
                    print(str(self.var_symbols[i]), np.std(x[i, :]))
            else:
                self.SteadyValues()
                try:
                    for i in range(self.n):
                        hp_cycle, hp_trend = sm.tsa.filters.hpfilter(100*np.log(x[i, :]+self.steady_values[self.steadys[i]]),
                                                                     lamb=hp_filter)
                        print(str(self.var_symbols[i]), '\t\t {:.4f}'.format(np.std(hp_cycle)))
                except:
                    print('Error: hp_filter takes only numbers as parameters')

            return None

        def stoch_simul(self, irf=40):

            self.TimeIteration()

            FT = self.F.T
            QT = self.Q.T

            print('\n'+'POLICY AND TRANSITION FUNCTIONS'+'\n')

            header = '\t'
            for v in self.var_symbols:
                header += '\t' + str(v)
            print(header)

            line = ''
            for i in range(self.n):
                line += '\t%.4f' % self.steady_values[self.steadys[i]]
            print('Constant' + line)

            for i in range(self.n):
                if (FT[i] != np.zeros((1, self.n))).any():
                    line = '\t'
                    for j in range(self.n):
                        line += '\t%.4f' % FT[i, j]
                    print(str(self.var_symbols[i])+'(-1)', line)
            for i in range(self.p):
                line = '\t'
                for j in range(self.n):
                    line += '\t%.4f' % QT[i, j]
                print(str(self.varexo_symbols[i]), '   ', line)

            # Impulse response functions
            if irf > 0:
                x = np.zeros((self.n, irf+2))
                ɛ = np.zeros((self.p, irf+2))
                ɛ[:, 1] = 1

                for t in range(1, irf+2):
                    x[:, t] = self.F @ x[:, t-1] + self.Q @ ɛ[:, t]

                for i in range(self.n):
                    plt.plot(x[i, 1:].T, 'k')
                    plt.hlines(0, 0, irf, 'r')
                    plt.title(str(self.var_symbols[i]))
                    plt.show()

            return None

In [30]:

    var = 'y c i k h w r z yh R'
    varexo = 'e'

    param_values = {'alpha':0.33, sy.symbols('beta'):0.99, 'delta':0.025, 'phi':1.75, 'sigma':1, 'rho':0.962159059423142}

    model = ('-log(z) + rho*log(z(-1)) + e',
             '-k + i + (1-delta)*k(-1)',
             '-c^(-sigma) + betta*c(+1)**(-sigma)*(1+r(+1))',
             '-y + z*k(-1)^alpha*h^(1-alpha)',
             '-r + alpha*y/k(-1) - delta',
             '-w + (1-alpha)*y/h',
             '-phi/(1-h) + w/c',
             '-y + c + i',
             '-yh + y/h',
             '-R + 1+r')

    initval = (1, 0.8, 0.2, 10, 0.33, 2, 0.01, 1, 3, 1.01)

In [31]:

    rbc = Dynare(var, varexo, param_values, model, initval)

    rbc.system

Out[31]:
$$\left[\begin{matrix}e_{t} + \rho \log{\left (z_{t-1} \right )} - \log{\left (z_{t} \right )}\\i_{t} + k_{t-1} \left(- \delta + 1\right) - k_{t}\\\beta c_{t+1}^{- \sigma} \left(r_{t+1} + 1\right) - c_{t}^{- \sigma}\\h_{t}^{- \alpha + 1} k_{t-1}^{\alpha} z_{t} - y_{t}\\\frac{\alpha y_{t}}{k_{t-1}} - \delta - r_{t}\\- w_{t} + \frac{y_{t}}{h_{t}} \left(- \alpha + 1\right)\\- \frac{\phi}{- h_{t} + 1} + \frac{w_{t}}{c_{t}}\\c_{t} + i_{t} - y_{t}\\- yh_{t} + \frac{y_{t}}{h_{t}}\\- R_{t} + r_{t} + 1\end{matrix}\right]$$
In [32]:

    rbc.steady()

    rbc.stoch_simul(irf=40)

    STEADY-STATE RESULTS

    y 	1.0058
    c 	0.7694
    i 	0.2364
    k 	9.4556
    h 	0.3336
    w 	2.0203
    r 	0.0101
    z 	1.0000
    yh 	3.0153
    R 	1.0101

    POLICY AND TRANSITION FUNCTIONS

    		y	c	i	k	h	w	r	z	yh	R
    Constant	1.0058	0.7694	0.2364	9.4556	0.3336	2.0203	0.0101	1.0000	3.0153	1.0101
    k(-1) 		0.0173	0.0437	-0.0264	0.9486	-0.0088	0.0881	-0.0031	0.0000	0.1315	-0.0031
    z(-1) 		1.3903	0.3396	1.0508	1.0508	0.2092	1.5257	0.0485	0.9622	2.2772	0.0485
    e     		1.4450	0.3529	1.0921	1.0921	0.2174	1.5857	0.0504	1.0000	2.3668	0.0504

In [33]:

    rbc.SimulatedMoments(hp_filter=1600, shocks_stderr=0.007828026718977966)

    SIMULATED MOMENTS

    VARIABLE 	 STD. DEV.
    y 		 1.4828
    c 		 0.5294
    i 		 4.8257
    k 		 0.4188
    h 		 0.6744
    w 		 0.8312
    r 		 5.3487
    z 		 1.0279
    yh 		 0.8312
    R 		 0.0523

In [34]:

    print('Standard Deviations')
    print(hp_cycles.std())

    Standard Deviations
    Output           1.602339
    Consumption      0.858648
    Investment       4.539281
    Capital          0.571145
    Hours            1.599009
    Wages            0.834119
    Interest Rate    0.392727
    TFP              0.996028
    Productivity     1.302317
    Price Level      0.888267
    dtype: float64