In the last blog, we went through the economic intuitions behind this consumption model, conducted some preliminary analysis in search for structural breaks and checked the stationarity of the endogenous variables using ADF and KPSS tests. Since we have found that none of the variables can be considered stationary, we can therefore ask whether the variables form a cointegrated system with a given number of “common trends”. Intuitively, I would argue that there exists one cointegration equation among the three variables as in the long run, one’s consumption and net wealth should add up to the one’s income.

So let’s keep going down to the path, and see how we can test for cointegration and include it in a VECM.

Cointegration Test

Given that all variables are consistent with the unit root hypothesis and possibly cointegrated, we can use Johansen’s trace statistic to test for the existence of a common trend i.e., a long-run cointegration relationship among the three endogenous variables.

First, we should estimate the VAR in levels so as to find out the optimal lags to use for the cointegration test. You might want to do this in first differences depending on your preferences, but basically the objective here is to run a preliminary VAR without worrying about whether the variables are cointegrated, which is obviously an issue given that they all have unit roots.

VAR in Levels

I’m going to allow for 10 lags in the first instance just to get the process going, and these are the variables in the VAR as I include the dummy variables sb_1975_4, sb_2008_3 to allow for the possibility of a structural breaks.

The BIC and HQ statistics suggest two lags for the optimal number of lags, whereas the AIC and FPE statistics suggest three. We choose the VAR model with three lags to allow for the possibility.

1

train = train.iloc[1:]

1
2
3

model = VAR(endog=train[['ln_rc','ln_rdy','ln_rnw']], exog =train[['sb_1975_4','sb_2008_3','const']])
res = model.select_order(10)
res.summary()

Johansen Cointegration Test

Then, we conduct the Johansen’s test with two lags. One has to ensure that the lag interval for the
differenced endogenous variables is 1 smaller than the number of lags used in VAR in levels by moving to the corresponding VECM form.

1
2
3

vec_rank = vecm.select_coint_rank(endog=train[['ln_rc','ln_rdy','ln_rnw']], det_order = 0, k_ar_diff = 2, method = 'trace', signif=0.05)
print(vec_rank.summary())
print('The rank to choose according to the Johansen test: '+ str(vec_rank.rank))

Johansen cointegration test using trace test statistic with 5% significance level
=====================================
r_0 r_1 test statistic critical value
-------------------------------------
  0   3          33.75          29.80
  1   3          11.44          15.49
-------------------------------------
The rank to choose according to the Johansen test: 1

The results suggest that we can reject the null hypothesis of no cointegration (or zero cointegrating vectors) using a 5% significance level. Moreover, we cannot reject the null hypotheses of at most 1 cointegrating vectors versus the alternative of 2. Therefore, we assume that there exists one (and only one) cointegrating vector.

VECM Estimation

Since there is a cointegration relationship existing, we should construct an ECM equation (e.g., an error correction model) that can explain the actual behavior of consumption during the sample period based on the long-run specification.

In developing the ECM model, we can include additional exogenous variables in the VECM to explain short-run dynamics of consumption. Besides the structural break around 1976 and the inequality gap, here are a bunch of other variables that we should consider. Remember the the determinants we discussed at the beginning?

• Disposable Income (in the current period) (+)
• Expectations (consumer confidence indexes, employment growth) (+)
• Wealth (stock market performance, housing prices) (+)
• Uncertainty (precautionary saving) (-)
• Availability of Credit (+)
• Real Interest Rate (?)
  • Substitution effect (-)
  • Income effect (+)

Explore Exogenous Variables

In developing the ECM model, we should include additional exogenous variables in the VECM to explain short-run dynamics of consumption.

Expectations

Here we use the unemployment rate and the consumer confidence index as exogenous variables. These variables are proxies for changes in the level of income uncertainty facing the household sector.

1
2
3
4
5

f, (ax1, ax2) = plt.subplots(1, 2,figsize=(14,4))
ax1.plot(train['unemp'])
ax2.plot(train['unemp'].diff().dropna())
ax1.set_title('unemp')
ax2.set_title('diff_unemp')

1
2
3
4
5

f, (ax1, ax2) = plt.subplots(1, 2,figsize=(14,4))
ax1.plot(train['consumer_confidence'])
ax2.plot(np.log(train['consumer_confidence']))
ax1.set_title('consumer_confidence')
ax2.set_title('ln_consumer_confidence')

Let’s then conduct the stationary tests for the change in the unemployment rate and the logarithm of the consumer confidence index, to check whether they are separately consistent with the stationarity assumption.

1
2

augmented_dickey_fuller_test(train['unemp'].diff().dropna())
kpss_test(train['unemp'].diff().dropna())

ADF Statistic: -5.032319
ADF p-value: 0.000019
stationary - null hypothesis of a unit root can be rejected at a 5% significant level
KPSS Statistic: 0.048204
KPSS p-value: 0.100000
stationary - null hypothesis of stationary cannot be rejected at a 5% significant level

1
2

augmented_dickey_fuller_test(np.log(train['consumer_confidence']))
kpss_test(np.log(train['consumer_confidence']))

ADF Statistic: -3.450606
ADF p-value: 0.009350
stationary - null hypothesis of a unit root can be rejected at a 5% significant level
KPSS Statistic: 0.134007
KPSS p-value: 0.100000
stationary - null hypothesis of stationary cannot be rejected at a 5% significant level

Availability of Credit

Evidently, the availability and cost of credit will also affect decisions of households to consume. So the behavior of monetary aggregates in terms of the availability of credit and the interest rates in terms of the cost of credit will matter for our projections for private consumption.

(Honestly, it would be ridiculous to ignore the crazy amount of money that the fed has been pumping into the economy since 2008 and its consequential impact. The hopeless belief in the Keynesian’s monetary policy.)

1
2
3
4
5

f, (ax1, ax2) = plt.subplots(1, 2,figsize=(14,4))
ax1.plot(train['m2'])
ax2.plot(np.log(train['m2']).diff())
ax1.set_title('real m2')
ax2.set_title('diff real m2')

1
2
3
4
5

f, (ax1, ax2) = plt.subplots(1, 2,figsize=(14,4))
ax1.plot(train['dff'])
ax2.plot(train['dff'].diff())
ax1.set_title('Federal Funds Effective Rate')
ax2.set_title('Diff Federal Funds Effective Rate')

1
2

augmented_dickey_fuller_test(np.log(train['m2']).diff().dropna())
kpss_test(np.log(train['m2']).diff().dropna())

ADF Statistic: -5.040178
ADF p-value: 0.000018
stationary - null hypothesis of a unit root can be rejected at a 5% significant level
KPSS Statistic: 0.172896
KPSS p-value: 0.100000
stationary - null hypothesis of stationary cannot be rejected at a 5% significant level

Let’s adding these new variables to the dataframe.

1
2
3
4
5
6
7

df['d_unemp'] = df['unemp'].diff()
df['ln_consumer_confidence'] = np.log(df['consumer_confidence'])
df['d_ln_m2'] = np.log(df['m2']).diff()
df['d_dff'] = df['dff'].diff()
df['spread_10y3m'] = df['yield_10y']-df['yield_3m']

train,test = df[1:cut], df[cut:]

Model Estimation

Based on all previous analysis, we can eventually run a VECM. Based on the long-run equation, we specify the form of the model with three endogenous variables log(rc), log(rdy) and log(rnw), and allow for five additional exogenous variables.

The lag intervals for the differenced endogenous variables are 2 given that we estimated the VAR with 3 lags in levels.

1
2

endog = ['ln_rc','ln_rdy','ln_rnw']
exog = ['sb_1975_4','sb_2008_3','d_unemp','ln_consumer_confidence', 'd_ln_m2', 'd_dff']

1
2
3
4

vec = vecm.VECM(endog = train[endog],
                exog =train[exog],
                k_ar_diff = 2, coint_rank = 1, deterministic = "co", missing= 'drop')
vecm_fit = vec.fit()

The coefficients in the long-run equation are both statistically significantly different from zero (see the first three tables below), and most importantly, consistent with our economic priors i.e., income and wealth have positive coefficients on consumption (see parameters for ln_rdy and ln_rnw in the first table).

For the exogenous variables, the negative coefficient on the structural breaks indicates that the increasing inequality may compress the consumption of households. Positive and significant coefficient on real m2 complies with our assumption that expanding monetary policy through quantitative ease may stimulate consumption in the short term.

Also, the coefficient on the ECM term is negative and significant with respect to real consumption, and significant and positive with respect to disposable income, both using a 5% level of significance, which are reasonable since both consumption and income are endogenous variables and dependent on each other, a rational household would react to disequilibrium in consumption relative to its long-run path though simultaneous adjustments - lower their spending (negative coefficient) and increase their income (positive coefficient).

From the “cointegration relations” table, we can verify our idea that in the long run, amount that one’s consumption and wealth accumulated would be equal to the income one earns.

1

vecm_fit.summary().as_latex

<bound method Summary.as_latex of <class 'statsmodels.iolib.summary.Summary'>
"""
Det. terms outside the coint. relation & lagged endog. parameters for equation ln_rc
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.6625      0.154     -4.310      0.000      -0.964      -0.361
exog1         -0.0027      0.001     -3.309      0.001      -0.004      -0.001
exog2         -0.0052      0.001     -4.940      0.000      -0.007      -0.003
exog3         -0.0069      0.001     -4.782      0.000      -0.010      -0.004
exog4          0.1386      0.032      4.279      0.000       0.075       0.202
exog5          0.1697      0.033      5.189      0.000       0.106       0.234
exog6          0.0006      0.000      1.463      0.144      -0.000       0.001
L1.ln_rc      -0.2099      0.068     -3.108      0.002      -0.342      -0.078
L1.ln_rdy      0.0742      0.043      1.739      0.082      -0.009       0.158
L1.ln_rnw      0.0485      0.022      2.248      0.025       0.006       0.091
L2.ln_rc      -0.0623      0.063     -0.986      0.324      -0.186       0.062
L2.ln_rdy      0.0077      0.042      0.184      0.854      -0.074       0.090
L2.ln_rnw      0.0039      0.022      0.176      0.860      -0.039       0.047
Det. terms outside the coint. relation & lagged endog. parameters for equation ln_rdy
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.2331      0.251     -0.928      0.353      -0.725       0.259
exog1         -0.0025      0.001     -1.879      0.060      -0.005       0.000
exog2         -0.0031      0.002     -1.777      0.076      -0.006       0.000
exog3         -0.0069      0.002     -2.937      0.003      -0.011      -0.002
exog4          0.0640      0.053      1.210      0.226      -0.040       0.168
exog5          0.1607      0.053      3.007      0.003       0.056       0.265
exog6         -0.0003      0.001     -0.445      0.656      -0.002       0.001
L1.ln_rc       0.1241      0.110      1.124      0.261      -0.092       0.340
L1.ln_rdy     -0.1826      0.070     -2.620      0.009      -0.319      -0.046
L1.ln_rnw      0.0519      0.035      1.469      0.142      -0.017       0.121
L2.ln_rc      -0.1472      0.103     -1.426      0.154      -0.350       0.055
L2.ln_rdy      0.0505      0.068      0.738      0.460      -0.084       0.184
L2.ln_rnw     -0.0868      0.036     -2.421      0.015      -0.157      -0.017
Det. terms outside the coint. relation & lagged endog. parameters for equation ln_rnw
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -1.4968      0.494     -3.028      0.002      -2.466      -0.528
exog1          0.0003      0.003      0.097      0.923      -0.005       0.005
exog2         -0.0009      0.003     -0.260      0.795      -0.007       0.006
exog3         -0.0024      0.005     -0.514      0.607      -0.011       0.007
exog4          0.3093      0.104      2.968      0.003       0.105       0.513
exog5          0.1839      0.105      1.748      0.080      -0.022       0.390
exog6         -0.0021      0.001     -1.576      0.115      -0.005       0.001
L1.ln_rc       0.1244      0.217      0.573      0.567      -0.301       0.550
L1.ln_rdy     -0.2364      0.137     -1.723      0.085      -0.505       0.033
L1.ln_rnw      0.0856      0.069      1.232      0.218      -0.051       0.222
L2.ln_rc       0.2163      0.203      1.065      0.287      -0.182       0.615
L2.ln_rdy     -0.2192      0.135     -1.629      0.103      -0.483       0.044
L2.ln_rnw      0.0711      0.071      1.008      0.314      -0.067       0.209
               Loading coefficients (alpha) for equation ln_rc                
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ec1           -0.0605      0.020     -2.993      0.003      -0.100      -0.021
               Loading coefficients (alpha) for equation ln_rdy              
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ec1            0.0887      0.033      2.684      0.007       0.024       0.153
               Loading coefficients (alpha) for equation ln_rnw              
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ec1           -0.1382      0.065     -2.126      0.033      -0.266      -0.011
          Cointegration relations for loading-coefficients-column 1          
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
beta.1         1.0000          0          0      0.000       1.000       1.000
beta.2        -1.1175      0.073    -15.311      0.000      -1.261      -0.974
beta.3         0.0555      0.058      0.963      0.336      -0.057       0.168
==============================================================================
""">

VECM Forecast

Out-of-sample Forecast

As we are comfortable with the model estimated, let’s do some forecasting with the out-of-sample data and evaluate how the model works.

To forecast for time beyond the sample period pre-2017, we need to provide futures shocks post-2017 i.e., the exogenous variables, which includes the dramatic economic downturn during the Covid recession. Let’s plot them out and get a feel.

1
2
3
4
5
6
7
8

f, (ax1, ax2, ax3) = plt.subplots(1,3,figsize=(18,4))
ax1.plot(test['d_unemp'])
ax2.plot(test['d_ln_m2'])
ax3.plot(test['d_dff'])
ax1.set_title('Diff Unemployment Rate')
ax2.set_title('Diff Real M2')
ax3.set_title('Diff Fed Fund Rate')
f.autofmt_xdate()

When forecasting, I include the confidence intervals for the predictions as well. Remember, point predictions won’t help you much with gauging the uncertainty of forecasts.

1
2
3
4
5
6
7
8
9

mid, lower_1se, upper_1se = vecm_fit.predict(steps=len(test),alpha = 0.32, exog_fc =test[exog])
_, lower_2se, upper_2se = vecm_fit.predict(steps=len(test),alpha = 0.05, exog_fc =test[exog])

pred_mid = pd.DataFrame(np.exp(mid), columns=['pred_rc','pred_rdy','pred_rnw'], index = test.index)
pred_lower_1se = pd.DataFrame(np.exp(lower_1se), columns=['lower_1se_rc','lower_1se_rdy','lower_1se_rnw'], index = test.index)
pred_upper_1se = pd.DataFrame(np.exp(upper_1se), columns=['upper_1se_rc','upper_1se_rdy','upper_1se_rnw'], index = test.index)
pred_lower_2se = pd.DataFrame(np.exp(lower_2se), columns=['lower_2se_rc','lower_2se_rdy','lower_2se_rnw'], index = test.index)
pred_upper_2se = pd.DataFrame(np.exp(upper_2se), columns=['upper_2se_rc','upper_2se_rdy','upper_2se_rnw'], index = test.index)
ptest = pd.concat([test, pred_mid, pred_lower_1se, pred_upper_1se, pred_lower_2se, pred_upper_2se],axis = 1)

The saving rate (percentage points) can be derived from total nominal consumption and
nominal household disposable income using the following formula:

$saving\ rate=100*(rdy - (rc + (government\ transfers + interest\ payments/gdp\ deflator))) /rdy$

1

ptest['pred_saving_rate']=100*(ptest.pred_rdy-(ptest.pred_rc+(ptest.gov_transfers+ptest.interest_payments/ptest.c_deflator)))/ptest.pred_rdy

Now let’s see how the model predicts,

1
2
3
4
5
6
7
8
9
10
11
12
13

def plot_point_prediction(test):
    fig, axes = plt.subplots(2, 2, figsize=(12,8))
    count=0
    for endog in ['rc','rdy','rnw','saving_rate']:
        col=count %2
        row=count //2
        ax=axes[row][col]
        ax.set_title(endog)
        ax.plot(test[endog], color='Black',linewidth=1, label="Observed")
        ax.plot(test['pred_'+endog],color='grey',linestyle='--',linewidth=1, label="Forecast")  
        ax.legend()
        count+=1
    fig.autofmt_xdate()

1

plot_point_prediction(ptest)

The results are not so bad. Actually pretty dope pre-pandemic. For a model that is trying to predict for consumption three years out, its predictive accuracy is out-of-expectation.

Though not surprisingly, the forecasts are over predicting actual real disposable income, real net wealth and therefore real consumption expenditures. This reflects the I(1) nature of the variables involved, the limited serial correlation in their innovation sequences, and the fact that all variables were trending upwards before the pandemic. Because of their I(1) nature, a dynamic forecast will effectively use the last known value (four years ago) of these variables to generate out-of-sample forecasts. Consequently, the forecasting model does not do good job of predicting the impact of covid19 on real consumption, and its predictions for the saving ratio appear to be systematically below the actual saving outcome.

1
2
3
4
5
6
7
8
9
10
11
12
13

def plot_forecast_with_confidence_intervals(df,test):
    fig, axes = plt.subplots(3, 1, figsize=(10,10))
    count=0
    for endog in ['rc','rdy','rnw']:
        col=count
        ax=axes[col]
        ax.set_title(endog)
        ax.plot(df[endog], color='Black',linewidth=1.2,label="Observed")
        ax.plot(test['pred_'+endog],'k--',linewidth=0.7,label="Forecast")
        ax.fill_between(test['pred_'+endog].index,test['lower_1se_'+endog],test['upper_1se_'+endog],alpha=0.07,color="k")
        ax.fill_between(test['pred_'+endog].index,test['lower_2se_'+endog],test['upper_2se_'+endog],alpha=0.07,color="k")      
        ax.legend(loc ='upper left')
        count+=1

With that being said, the previous assessment considers only a single point prediction using the baseline model. A fairer assessment of the model can be made considering the confidential intervals for prediction, which suggest that the underlying model is useful for predicting the outcome of the consumption and net wealth.

1

plot_forecast_with_confidence_intervals(df.iloc[200:],ptest)

While there is obvious (downward) bias in the model’s predictions, it should be noted that this partly a reflection of the initial starting point of the forecast. If one begins the forecast after the pandemic, the resulting confidential intervals would be much closer to what observed. In this way, a static forecast which is one-period ahead forecasting would do a much better job. This is also how the model should be used in practice. As always, predicting future is difficult, trying to predict mutiple years out is meaningless.

Long-run Forecasts Based on Potential Paths of Monetary Policy

Now it has come to a very interesting point. Let’s fit a VECM with the complete data from 1962 to 2021Q2, and use it to forecast for 2021Q3 onwards using two simulated forward paths of the Fed’s monetary policy. Correct assumptions would help us to get more accurate estimates for the long-run consumption, disposal income and net wealth of the American households after the pandemic.

1
2
3
4
5
6

endog = ['ln_rc','ln_rdy','ln_rnw']
exog = ['sb_1975_4','sb_2008_3','d_unemp','ln_consumer_confidence', 'd_ln_m2','d_dff']
vec = vecm.VECM(endog = df[endog],
                exog =df[exog],
                k_ar_diff = 2, coint_rank = 1, deterministic = "co", missing= 'drop')
vecm_fit = vec.fit()

Let’s generate two paths for the forward monetary shocks. One path assumes that fed keeps its current pacing of assets purchasing i.e., no tapering, the other path assumes four interest rate hikes in 2022, 2023 and 2024 assuming the same pace as the period after the financial crisis from 2014Q4 to 2017Q3, with all the other exogenous variables being the same, specifically, current value for unemployment rate and the historical average for consumer confidence.

1
2
3
4
5
6
7
8
9
10
11
12
13
14

forecast_range = pd.date_range(start='04/01/2021', end ='01/01/2024',freq='QS' )
shocks = {}
shocks = pd.DataFrame(shocks, index = forecast_range, columns = exog)
shocks['sb_1975_4'] = 1
shocks['sb_2008_3'] = 1
shocks['d_unemp'] = 0
shocks['ln_consumer_confidence'] = df['ln_consumer_confidence'].mean()

no_hikes_shocks = shocks.copy()
no_hikes_shocks['d_ln_m2'] = df.iloc[-4:]['d_ln_m2'].mean()
no_hikes_shocks['d_dff'] = 0.0

four_hikes_shocks =shocks.copy()
four_hikes_shocks[['d_ln_m2','d_dff']] = df.loc['2014-10-01':'2017-07-01',['d_ln_m2','d_dff']].values

Here’s how the money aggregates is going to be like if assuming no tapering until 2024, which would be of course totally insane…

1
2

# change of log(real m2) and change of fed fund rate
no_hikes_shocks[['d_ln_m2','d_dff']]

Here’s a more likely path assumes that the fed would stop expanding its balance sheet from the forth quarter of 2021 and hike once in 2022 and three times at the end 2023 and beginning of 2024. This path is quite aggressive.

(By the way, when central banks expend their balance sheets i.e., increase their net asset position, they literally buy assets or reduce their liability through money creation, which means writing checks to themselves.)

1
2
3

# hike schedule - historical simulation using 2014-2017 data
## change of log(real m2) and change of fed fund rate
four_hikes_shocks[['d_ln_m2','d_dff']]

1
2
3
4
5
6
7
8
9
10
11
12

mid_no_hikes = vecm_fit.predict(steps=12,exog_fc =no_hikes_shocks[exog])
pred_mid_no_hikes = pd.DataFrame(np.exp(mid_no_hikes), columns=['pred_rc','pred_rdy','pred_rnw'], index = forecast_range)

mid, lower_1se, upper_1se = vecm_fit.predict(steps=12,alpha = 0.32, exog_fc =four_hikes_shocks[exog])
_, lower_2se, upper_2se = vecm_fit.predict(steps=12,alpha = 0.05, exog_fc =four_hikes_shocks[exog])

pred_mid = pd.DataFrame(np.exp(mid), columns=['pred_rc','pred_rdy','pred_rnw'], index = forecast_range)
pred_lower_1se = pd.DataFrame(np.exp(lower_1se), columns=['lower_1se_rc','lower_1se_rdy','lower_1se_rnw'], index = forecast_range)
pred_upper_1se = pd.DataFrame(np.exp(upper_1se), columns=['upper_1se_rc','upper_1se_rdy','upper_1se_rnw'], index = forecast_range)
pred_lower_2se = pd.DataFrame(np.exp(lower_2se), columns=['lower_2se_rc','lower_2se_rdy','lower_2se_rnw'], index = forecast_range)
pred_upper_2se = pd.DataFrame(np.exp(upper_2se), columns=['upper_2se_rc','upper_2se_rdy','upper_2se_rnw'], index = forecast_range)
pred = pd.concat([pred_mid, pred_lower_1se, pred_upper_1se, pred_lower_2se, pred_upper_2se],axis = 1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

def plot_forecasts_given_different_paths(df,pred1,pred2):
    fig, axes = plt.subplots(3, 1, figsize=(10,10))
    count=0
    for endog in ['rc','rdy','rnw']:
        col=count
        ax=axes[col]
        ax.set_title(endog)
        ax.plot(df[endog], color='Black',linewidth=1.2,label="Observed")
        ax.plot(pred1['pred_'+endog],'--',color='tab:blue', linewidth=1.5,label="Forecast - no hikes")
        ax.plot(pred2['pred_'+endog],'--',color='tab:red', linewidth=1.5,label="Forecast - four hikes")
        ax.fill_between(pred2['pred_'+endog].index,pred2['lower_1se_'+endog],pred2['upper_1se_'+endog],alpha=0.07,color="tab:red")
        ax.fill_between(pred2['pred_'+endog].index,pred2['lower_2se_'+endog],pred2['upper_2se_'+endog],alpha=0.07,color="tab:red")
        ax.legend(loc ='upper left')
        ax.axvline(x=datetime(2022, 7, 1), color='tab:red')
        ax.axvline(x=datetime(2023, 7, 1), color='tab:red')
        ax.axvline(x=datetime(2023, 10, 1), color='tab:red')
        ax.axvline(x=datetime(2024, 1, 1), color='tab:red')
        count+=1

In March 2020, in response to the economic impact of the pandemic, the Fed acted quickly having learned its lesson from the financial crisis and cut short-term interest rates to zero. It also restarted its large-scale asset purchases and since June 2020, the Fed has been buying 80 billion of Treasury securities and 40 billion of mortgage-backed securities each month, with its balance sheet swelling from 4.4 trillion to 8.6 trillion. The extreme monetary support combined with a huge fiscal stimulus package holds the U.S economy and shields American household incomes. Though plots won’t be shown here, the model forecasts a much lower rdy and rnw if not taking into account the rapid increase of money base. Therefore, in the last 20 months, we all see the quick recovery of consumption and a rapid growth of asset prices, which is also reflected in the steeper slope of the third graph. However, as the economy rebounded, it may no longer need such extreme measures of support and keeping them in place could do more harm than good. For example, low mortgage rates have fueled a boom in house prices, but the problems that now afflict the economy are mostly supply issues while demand, which the bond buys most directly affects. Therefore, the policy makers may start to hit the brakes. That’s also why we are seeing Fed officials talking about tapering the pace of its bond purchases and lately laying the groundwork for higher rates.

1. Growth rates return to normal

First it’s not surprising that real consumption, disposable income and household net wealth are projected to grow at milder rates comparing to those when the economy was bottoming out from the slump. As the stimulus aids wear off over time, one would expect the growth rates to normalize to pre-pandemic levels.

2. Growth differences taking expansionary and tight monetary policy

Secondly, assuming a tighter monetary policy to come, we project that the growth rate of real consumption throughout 2023 would slow down to 2.2% from 3.0%. The real income and net wealth of households would decrease from 3.9% and 9.1% to 2.4% and 6.5% respectively. At the end of 2023, all three measures would be lower than those come from the first path - that the Fed sticks to its current rate target and asset buying pace. In fact the amounts would be lowered by 1.8%, 3.3% and 4.8% respectively.

3. Foreseeable shallow decline in mid 2022

Thirdly, we also notice that in the middle of 2022, the private consumption and household income would show a temporary and shallow decline and then grow at an annual rate again coming close to the path observed just before the pandemic. Likely this is caused by the tightening of policy in 2022. So far, Powell has announced that the fed will begin tapering this month, the first step toward pulling back on the massive amount of help it had been providing markets and the economy. One the fiscal side, though public programs to support households and businesses have been the most powerful engine of recovery from the hit, governments of the major economies now are hitting the brakes. The money they’ ll pull out of their economies in 2022 amounts to 2.5 percentage points of the world’s GDP, five times bigger than anything that happened during the turn to austerity after the 2008 crisis, according to UBS estimates. The fiscal tightening - the withdrawal of pandemic spending will likely have more impact on the global economy next year.

1

plot_forecasts_given_different_paths(df.iloc[220:],pred_mid_no_hikes,pred)

It’s official on Nov.3rd that the Federal Reserve is winding down its aggressive pandemic-era asset purchasing, or as what we call “tapering”. The announcement actually didn’t come out as a surprise since Fed officials had been laying the groundwork for policy tightening as the economy gradually recovered from the Covid19 slump. Furthermore, due to the lately surging inflation readings, we are seeing traders start to bet on faster Fed hikes. So it seems interesting to get some ideas of how the post-pandemic consumption behaviors of the American households would be, and quantify the potential impact of policy shifting on these behaviors using two classic time series models - the Vector Autoregression (VAR) and the Vector Error Correction Model (VECM). Tribute to Christopher Sims who first promoted the use of VAR in empirical macroeconomics.

We will build the models following the steps below,

And before we get onto that, I will introduce some of the economic theories behind predicting real consumption.

Economic Intuitions Behind

Let’s first recap some of the main determinants of private consumption,

• Disposable Income (in the current period) (+)
• Expectations (consumer confidence indexes, employment growth) (+)
• Wealth (stock market performance, housing prices) (+)
• Uncertainty (precautionary saving) (-)
• Availability of Credit (+)
• Real Interest Rate (?)
      - Substitution effect (-)
      - Income effect (+)

You might remember that perhaps the simplest theories of consumption, Keynesian theories of consumption, would link private consumption to disposable income in the current period, and of course positively. So as disposable income goes up, private consumption is expected to go up as well. More sophisticated, forward looking models would give more prominence or link consumption not to disposable income today but to expectations about the evolution of disposable income in the future. A number of indicators that can be used to gauge these expectations, for example, consumer confidence index, employment growth or my personally preferred - yield spreads. These can be used as a proxy for expectations for future income.

Another variable thought to be important by economists in terms of explaining the behavior of private consumption is wealth. In that sense, for this wealth channel, the performance of the stock market or even the behavior of housing prices, would have an impact on private consumption.

Another important factor is uncertainty, which would actually have a negative association with private consumption because it would create incentives for what the economists called precautionary saving. So agents will try to create a buffer to be able to fall back on in case a bad shock hits.

The availability of credit is also thought to be an important determinant with, obviously, a positive relationship.

And finally, the real interest rate is also thought to be one of the important determinants of private consumption.
The effect of real interest on consumption can be a little bit ambiguous. Perhaps the most straightforward or direct effect
would be what economists call a substitution effect where, as interest rates go up, the incentives to consume are reduced. So agents will have, on the contrary, incentives to save more. And therefore, an increase in rates will be linked to a decrease in private consumption. But on the other hand, it is also possible that if consumers are net savers the increase in interest rates
might generate some additional income– what economists call an income effect–and this might actually lead to an increase in private consumption.

Getting Data

Based on the economic intuitions, I grabbed the following well-defined and organized economic data from the Federal Reserve Economic Data of the St. Louis Fed. The FRED allows modelers to write programs and build applications that retrieve data through an API provided. There is a list of packages/libraries that one can turn to, depending on one’s programming language. Since I have to model everything in Python as it’s compatible with the github blog, I use the library fredapi to request the data. However, in practice, I would probably go for EViews as most economists prefer and do.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

import pandas as pd
import numpy as np
from fredapi import Fred
import matplotlib.pyplot as plt
from datetime import datetime
%matplotlib inline

import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.stattools import kpss
from statsmodels.tsa.api import VAR
from statsmodels.tsa.vector_ar import vecm
from statsmodels.tsa import filters

pd.options.display.max_colwidth = 50
pd.options.display.max_rows = 400

1

fred = Fred(api_key=your_key_string)

The period I study is from 1962Q1 to 2021Q2, which ensures that the economic data that I am interested in all has a value.

1
2

startDate = '1962-01-01'
endDate = '2021-04-01'

1
2
3
4

def get_series_and_print_info(fred, series_id, startDate, endDate, name):
    series = fred.get_series(series_id, observation_start=startDate, observation_end=endDate, frequency='q')
    print(name+': '+fred.get_series_info(series_id).title)
    return series

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

data = {}
data['saving_rate'] = get_series_and_print_info(fred, 'a072rc1q156sbea',startDate, endDate, 'saving_rate')
data['ns'] = get_series_and_print_info(fred, 'w986rc1q027sbea',startDate, endDate, 'ns')
data['nnw'] = get_series_and_print_info(fred, 'TNWBSHNO',startDate, endDate, 'nnw')
data['nl'] = get_series_and_print_info(fred, 'TCMILBSHNO',startDate, endDate, 'nl')
data['na'] = get_series_and_print_info(fred, 'TABSHNO',startDate, endDate, 'na')
data['nfa'] = get_series_and_print_info(fred, 'TFAABSHNO',startDate, endDate, 'nfa')
data['interest_payments'] = get_series_and_print_info(fred, 'b069rc1q027sbea',startDate, endDate, 'interest_payments')
data['cpi'] = get_series_and_print_info(fred, 'cpiaucsl',startDate, endDate, 'cpi')
data['c_deflator'] = get_series_and_print_info(fred, 'dpcerd3q086sbea',startDate, endDate, 'c_deflator')/100
data['dff'] = get_series_and_print_info(fred, 'DFF',startDate, endDate, 'dff')
data['yield_3m'] = get_series_and_print_info(fred, 'dgs3',startDate, endDate, 'yield_3m')
data['yield_10y'] = get_series_and_print_info(fred, 'irltlt01usq156n',startDate, endDate, 'yield_10y')
data['prime_rate'] = get_series_and_print_info(fred, 'dprime',startDate, endDate, 'prime_rate')
data['m2'] = get_series_and_print_info(fred, 'M2REAL',startDate, endDate, 'm2')
data['ndy'] = get_series_and_print_info(fred, 'dpi',startDate, endDate, 'ndy')
data['rdy'] = get_series_and_print_info(fred, 'dpic96',startDate, endDate, 'rdy')
data['npiy'] = get_series_and_print_info(fred, 'b703rc1q027sbea',startDate, endDate, 'npiy')
data['ny'] = get_series_and_print_info(fred, 'pincome',startDate, endDate, 'ny')
data['ry'] = get_series_and_print_info(fred, 'rpi',startDate, endDate, 'ry')
data['nc'] = get_series_and_print_info(fred, 'pcec',startDate, endDate, 'nc')
data['rc'] = get_series_and_print_info(fred, 'pcecc96',startDate, endDate, 'rc')
data['nc_services'] = get_series_and_print_info(fred, 'pcesv',startDate, endDate, 'nc_services')
data['nc_durable_goods'] = get_series_and_print_info(fred, 'pcdg',startDate, endDate, 'nc_durable_goods')
data['unemp'] = get_series_and_print_info(fred, 'unrate',startDate, endDate, 'unemp')
data['house_prices'] = get_series_and_print_info(fred, 'ussthpi',startDate, endDate, 'house_prices')
data['mortgage_30y'] = get_series_and_print_info(fred, 'MORTGAGE30US',startDate, endDate, 'mortgage_30y')
data['consumer_confidence'] = get_series_and_print_info(fred, 'CSCICP03USM665S',startDate, endDate, 'consumer_confidence')
data['gov_debt'] = get_series_and_print_info(fred, 'fygfdpun',startDate, endDate, 'gov_debt')
data['gov_transfers'] = get_series_and_print_info(fred, 'w211rc1q027sbea',startDate, endDate, 'gov_transfers')

    saving_rate: Personal saving as a percentage of disposable personal income
    ns: Net private saving: Households and institutions
    nnw: Households and Nonprofit Organizations; Net Worth, Level
    nl: Households and Nonprofit Organizations; Debt Securities and Loans; Liability, Level
    na: Households and Nonprofit Organizations; Total Assets, Level
    nfa: Households and Nonprofit Organizations; Total Financial Assets, Level
    interest_payments: Personal interest payments
    cpi: Consumer Price Index for All Urban Consumers: All Items in U.S. City Average
    c_deflator: Personal consumption expenditures (implicit price deflator)
    dff: Federal Funds Effective Rate
    yield_3m: Market Yield on U.S. Treasury Securities at 3-Year Constant Maturity
    yield_10y: Long-Term Government Bond Yields: 10-year: Main (Including Benchmark) for the United States
    prime_rate: Bank Prime Loan Rate
    m2: Real M2 Money Stock
    ndy: Disposable Personal Income
    rdy: Real Disposable Personal Income
    npiy: Personal income receipts on assets: Personal dividend income
    ny: Personal Income
    ry: Real Personal Income
    nc: Personal Consumption Expenditures
    rc: Real Personal Consumption Expenditures
    nc_services: Personal Consumption Expenditures: Services
    nc_durable_goods: Personal Consumption Expenditures: Durable Goods
    unemp: Unemployment Rate
    house_prices: All-Transactions House Price Index for the United States
    mortgage_30y: 30-Year Fixed Rate Mortgage Average in the United States
    consumer_confidence: Consumer Opinion Surveys: Confidence Indicators: Composite Indicators: OECD Indicator for the United States
    gov_debt: Federal Debt Held by the Public
    gov_transfers: Personal current transfer payments

1

df = pd.DataFrame(data)

Since we are interested in predicting the real consumption, which means eliminating the effect of prices, nominal variables including consumption, disposable income and net wealth need to be converted into real terms. Here they are deflated by the price deflator for consumption (c_deflator).

1
2
3
4

df['rnw'] = df['nnw']/df['c_deflator']
df['rl'] = df['nl']/df['c_deflator']
df['ra'] = df['na']/df['c_deflator']
df['rfa'] = df['nfa']/df['c_deflator']

Data Visualization and Preliminary Analysis (Some Fun Ideas)

Structural Breaks

To get a feel for the data, let’s plot log(rc/rdy) against log(rnw/rdy). Also regress log(rc/rdy) on a constant and log(rnw/rdy) and then examine the fitted residuals, looking for evidence of any structural changes in the relationship between consumption, disposable income and net wealth.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

fig, ax1 = plt.subplots()
ax1.set_xlabel('Year')
ax1.set_ylabel('rdy rl')
ax1.plot(df['rc'], color='tab:grey',label='real consumption')
ax1.plot(df['rdy'], color='tab:blue',label='real disposable income')
ax1.plot(df['rl'], color='tab:green',label='real liability')
ax1.legend(loc='upper left')
ax2 = ax1.twinx()  
ax2.set_ylabel('rnw ra rfa') 
ax2.plot(df['rnw'], color='tab:red',label='real net worth')
ax2.plot(df['ra'], color='tab:orange',label='real asset')
ax2.plot(df['rfa'], color='tab:purple',label='real fincl asset')
ax2.legend(loc='lower right')
fig.tight_layout() 
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

data1 = df['rc']/df['rdy']
data2 = df['rc']/df['rnw']
data3 = df['rc']/df['rl']
data4 = df['rc']/df['ra']
data5 = df['rc']/df['rfa']
fig, ax1 = plt.subplots()
ax1.set_xlabel('Year')
ax1.set_ylabel('rdy')
ax1.plot(data1, color='tab:blue',label='rc/rdy')
ax1.legend(loc='upper left')
ax2 = ax1.twinx()  
ax2.set_ylabel('rnw ra rfa') 
ax2.plot(data2, color='tab:red',label='rc/rnw')
ax2.plot(data4, color='tab:orange',label='rc/ra')
ax2.plot(data5, color='tab:purple',label='rc/rfa')
ax2.legend(loc='upper right')
fig.tight_layout() 
plt.show()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

data1 = np.log(df['rc']/df['rdy'])
data2 = np.log(df['rnw']/df['rdy'])
fig, ax1 = plt.subplots()
color = 'tab:green'
ax1.set_xlabel('Year')
ax1.set_ylabel('rc/rdy', color=color)
ax1.plot(data1, color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax2 = ax1.twinx()  
color = 'tab:blue'
ax2.set_ylabel('rnw/rdy', color=color) 
ax2.plot(data2, color=color)
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout()
ax1.axvline(x=datetime(1975, 10, 1), color='k')
ax1.axvline(x=datetime(2008, 7, 1), color='k')
plt.show()

Are there any structural breaks in the relationship before 1975Q4 and after 2008Q3? If so, we have to create dummies to allow for them. Let’s run a regression to further confirm it.

1
2
3
4
5

Y = np.log(df['rc']/df['rdy'])
X = np.log(df['rnw']/df['rdy'])
X = sm.add_constant(X)
model = sm.OLS(Y,X)
res = model.fit()

1
2
3
4

ax = res.resid[:-8].plot(style='.',grid=True)
ax.xaxis.grid(True, linestyle='--', linewidth=0.25)
ax.axvline(x=datetime(1975, 10, 1), color='k')
ax.axvline(x=datetime(2008, 7, 1), color='k')

Both the time plot and the regression of log(rc/rdy) on a constant and log(rnw/rdy) suggests that there were structural breaks around 1975Q4 and 2008Q3. Notice the consistently negative residuals till 1975Q4 and the wider gap between log(rc/rdy) and log(rnw/rdy) before 1976 compared to the gap after 1976. Another one is that household net worth and assets have increased at an accelerating rate after the 2008 global financial crisis (increasing rnw/rdy ratio). While at the same time, we don’t see a rising trend of consumption compared to income even though households have accumulated higher real net wealth benefited from the booming asset prices post-crisis, which is a quite different behavior compared to pre-crisis.

These are strong suggestive evidence that structural breaks occurred around 1975Q4 and 2008Q3. But why? What are the drivers?

Inequality and the Kondratiev theory

Let’s now take a look at the inequality data, which will give us better clues to our questions.

Data is downloaded from the World Inequality Database WID. Buckets I selected are top 10%, middle 40% and the bottom 50% of the US population.

1
2
3
4

income_inequality = pd.read_csv('income_inequality.csv',index_col=0)
wealth_inequality = pd.read_csv('wealth_inequality.csv',index_col=0)
pretax_income = pd.read_csv('pretax_income.csv',index_col=0)
net_wealth = pd.read_csv('net_wealth.csv',index_col=0)

As we can tell from the following graphs, income inequality in the US has risen dramatically since the mid 1970s, see how the income of the top 10% grows whereas that of the middle and bottom classes remains nearly flat. This issue has drawn heightened attention in recent years (including me :D).

In the past decade, economic observers have also become increasingly worried about “secular stagnation” - or a chronic shortfall of aggregate demand, fearing that this shortfall will constrain American economic growth in coming years. These two phenomena—rising inequality and chronic weakness of demand - are related. Specifically, since the marginal propensity to consume (MPC) is much lower at the higher wealth quintiles (for low-wealth households, the MPC is 10 times larger than it is for wealthy households, see this paper elaborating it), rising inequality transfers income from high MPC households in the bottom and middle of the income distribution to lower MPC households at the top. All else equal, this redistribution away from low- to high-saving households reduces consumption spending, which drags on demand growth.

1

pretax_income.iloc[:,1:].plot(title='Pre-tax Income')

1

net_wealth.iloc[:,1:].plot(title='Net Wealth')

    <matplotlib.axes._subplots.AxesSubplot at 0x23c76c39e48>

Due to the differences in income growth, we see a larger and larger wealth gap between groups.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

data1 = wealth_inequality['sw_p99p100']
data2 = wealth_inequality['sw_p90p100']
data3 = wealth_inequality['sw_p50p90']
data4 = wealth_inequality['sw_p0p50']

fig, ax1 = plt.subplots()
color = 'tab:blue'
ax1.set_xlabel('Year')
ax1.set_ylabel('Rich', color=color)
ax1.plot(data1, color='royalblue',label='Share of p99p100')
ax1.plot(data2, color='dodgerblue',label='Share of p90p100')
ax1.tick_params(axis='y', labelcolor=color)
ax1.legend(loc ='center left')
ax1.set_title('Wealth Inequality')
ax2 = ax1.twinx()  
color = 'tab:green'
ax2.set_ylabel('Poor', color=color) 
ax2.plot(data3, color='limegreen',label='Share of p50p90')
ax2.plot(data4, color='lightgreen',label='Share of p0p50')
ax2.tick_params(axis='y', labelcolor=color)
fig.tight_layout() 
ax2.legend(loc ='center right')
plt.show()

Based on the theory, we can assume that there is a negative correlation between income inequality and consumption. Here I created an inequality indicator using the difference in shares of wealth of the top 10% and the bottom 50%, so that we can further look into the change of inequality over the last half century. The long-term trend is estimated using the Hodrick–Prescott filter, which can smooth the data by getting rid of the short-term fluctuations (shorter cycles).

1
2
3
4
5

wealth_inequality.index = pd.to_datetime(wealth_inequality.index, format='%Y')
# create a variable nw_inequlity for the wealth gap between top 10% and bottom 50%
wealth_inequality['nw_inequality'] = wealth_inequality['sw_p90p100']-wealth_inequality['sw_p50p90']
# apply the Hodrick-Prescott filter to smooth the inequality data and estimate the long-term trend
wealth_inequality['f_nw_inequality'] = pd.DataFrame(filters.hp_filter.hpfilter(wealth_inequality['nw_inequality'],lamb=6.25)[1])

1
2
3
4
5

f, (ax1, ax2) = plt.subplots(1,2,figsize=(15,4))
ax1.plot(wealth_inequality['nw_inequality'])
ax2.plot(wealth_inequality['f_nw_inequality'])
ax1.set_title('nw inequality')
ax2.set_title('filtered nw inequality (trend)')

Seeing the long-term inequity trend, it’s no wonder why we have identified two structural breaks.

Actually up to now, if you are familiar with the theory of Kondratieff cycles, this probably starts to look quite interesting to you. Kondratiev waves (also called super-cycles, great surges, long waves, K-waves or the long economic cycle) are hypothesized cycle-like phenomena in the modern world economy. The phenomenon is closely connected with the technology life cycle. It is stated that the period of a wave ranges from forty to sixty years, the cycles consist of alternating intervals of high sectoral growth and intervals of relatively slow growth.

Though there is a lack of agreement about both the cause of the waves and the start and end years of particular waves, inequity appears to be the most obvious driver, and yet some researchers have presented a technological and credit cycle explanation as well. However, among several modern timing versions of the cycles based on any of the causes, the consensus is that the start of our current cycle - Wave of the Information and Telecommunications Revolution falls into the range of early 1970s to mid 1980s, which happens to be the first break point we identified.

Every wave of innovations lasts approximately until the profits from the new innovation or sector fall to the level of other, older, more traditional sectors. It is a situation when the new technology, which originally increased a capacity to utilize new sources from nature, reached its limits and it is not possible to overcome this limit without an application of another new technology. For the end of an application phase of any wave there is typically an economic crisis and economic stagnation. The financial crisis of 2007–2008 is viewed by some as a result of the coming end of the Wave of the Information and Telecommunications Revolution, our second break point.

So after all, I create two dummy variables below to represent the turning points of the cycle.

1
2

df['sb_1975_4'] = np.where(df.index < '1976-10-1',0,1)
df['sb_2008_3'] = np.where(df.index < '2008-07-1',0,1)

Long-run Model Specification

Now let’s put aside the fun ideas and start model development with specifying the long-run predictive equation for real consumption.

Consider the following long-run model for U.S. real consumption:

\begin{align}
\log rc_t &= \beta_0 + \beta_1 \log rdy_t + \beta_1 \log rnw_t + \epsilon_t
\end{align}

It’s helpful to express our model in logarithmic form, so that the parameter estimates are elasticities.

1
2
3
4

df['ln_rc'] = np.log(df['rc'])
df['ln_rdy'] = np.log(df['rdy'])
df['ln_rnw'] = np.log(df['rnw'])
df['const'] = 1

Sample and Forecast Period

The sample period for model specification i.e., specifying the preferred model for forecasting purposes can be 1962:1 to 2016:4 (and no later than 2007:4). Remember we have to stick with the train data from now on for the model development. The test set will be preserved for testing the predictive power of the model.

1
2
3

cut = 220
train,test = df[:cut], df[cut:]
train.shape,test.shape

    ((220, 39), (18, 39))

1

train.tail()

Check for Stationary

I’m going to pre-test the data to determine whether if it has a unit root, that is, whether it’s stationary or non-stationary.

1. What is stationarity?

A stationary process is a stochastic process whose probability distribution does not change over time. For our work in time series analysis, we mostly care about stationarity in its weak form - covariance stationary.

Covariance stationary:

• Constant mean:
  \begin{align}
  &E(y_t) = E(y_{t+1}) = \mu
  \end{align}
• Constant variance:
  \begin{align}
  &Var(y_t) = Var(y_{t+1}) = \sigma^2
  \end{align}
• Covariance depends on time that has elapsed between observations, not on reference period:
  \begin{align}
  &Cov(y_t,y_{t+j}) = Cov(y_s,y_{s+j}) = \gamma_j
  \end{align}

To get an intuition of stationarity, one can imagine a frictionless pendulum. It swings back and forth in an oscillatory motion, yet the amplitude and frequency remain constant. Although the pendulum is moving, the process is stationary as its “statistics” are constant (frequency and amplitude). However, if a force were to be applied to the pendulum, either the frequency or amplitude would change, thus making the process non-stationary.

2. Why is stationarity important?

Stationarity is one important assumption underlying in time series analysis as if the series is non-stationary, the shocks do not die out, the impact of disturbances become more influential over time and so gives a poor forecast accuracy.

Consequences of non-stationarity

• The effect of a shock will diminish as time elapses
• Statistical consequences
       - Non-normal distribution of test statistics
       - Bias in autoregressive coefficients; we might mistakenly estimate an AR(1), deficient forecast
       - Usual confidence intervals for coefficients not valid, poor forecast ability

3. How do we determine whether a time series is nonstationary?

As always, visual inspection is one of the most useful ways to identify non stationary series. It looks that all of the curves obviously have a trend, which is definitely against the definition of stationarity.

1

train[['ln_rc','ln_rdy','ln_rnw']].plot()

    <matplotlib.axes._subplots.AxesSubplot at 0x23c0007c248>

Then we can turn to statistical tests. The first would be Augmented Dickey-Fuller (ADF) Test, which tests for non-stationarity.

Recall AR(1) model:
\begin{align}
y_t =  by{t_1} + \epsilon_t
\end{align}

A special case is the random walk when $b = 1$. However, for stationarity, it requires $b<1$.

If generalizing to AR(p), it means that roots of the polynomial below must all be >1 in abs value,
\begin{align}
1 - b_1z -b_2z^2-b_3z^3-…-b_pz^p
\end{align}

If one of the roots = 1, then it is said to have a unit root i.e., non-stationary.

The ADF Test can test for the coefficient on $y_{t-1}$. Basically, it regresses y on its lag testing for significance of coefficient. I won’t go into great mathematical details, if anyone’s interested, check out this wiki page.

However, ADF has been found to have low power in certain circumstances:

• stationary processes with near-unit roots. For example, it has difficulty distinguishing between b = 1 and b = 0.95, especially with small samples.
• Trend stationary processes.

So alternative tests have been designed, which are Phillips–Perron (PP) Test and Kwiatkowski–Phillips–Schmidt–Shin (KPSS)Test. I usually prefer to use all three of them together, but here I would focus on ADF and KPSS as they have been implemented in statsmodels.

A few notes:
 - ADF and PP are called unit root tests; the null hypothesis is that $y_t$ has a unit root; is I(1) or higher.
 - KPSS, on the other hand, is a stationarity test; the null hypothesis is that $y_t$ is I(0).
 - Correct specification is key. Intercept and trend should be included when appropriate as critical values for the t-statistics will vary depending on whether they are included.
 - Structural breaks can complicate matters further. That’s why we start this blog with detecting structural breaks.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

def augmented_dickey_fuller_test(time_series):
    result = adfuller(time_series.values)
    print('ADF Statistic: %f' % result[0])
    print('ADF p-value: %f' % result[1])
    if result[1]<0.05:
        print('stationary - null hypothesis of a unit root can be rejected at a 5% significance level')
    else:
        print('non-stationary - null hypothesis of a unit root cannot be rejected at a 5% significance level')

def kpss_test(timeseries):
    result = kpss(timeseries, regression="c")
    print('KPSS Statistic: %f' % result[0])
    print('KPSS p-value: %f' % result[1])
    if result[1]>=0.05:
        print('stationary - null hypothesis of stationary cannot be rejected at a 5% significance level')
    else:
        print('non-stationary - null hypothesis of stationary can be rejected at a 5% significance level')

1
2

augmented_dickey_fuller_test(train['ln_rc'])
kpss_test(train['ln_rc'])

    ADF Statistic: -1.949803
    ADF p-value: 0.309007
    non-stationary - null hypothesis of a unit root cannot be rejected at a 5% significance level
    KPSS Statistic: 2.297233
    KPSS p-value: 0.010000
    non-stationary - null hypothesis of stationary can be rejected at a 5% significance level

1
2

augmented_dickey_fuller_test(train['ln_rdy'])
kpss_test(train['ln_rdy'])

    ADF Statistic: -2.970069
    ADF p-value: 0.037786
    stationary - null hypothesis of a unit root can be rejected at a 5% significance level
    KPSS Statistic: 2.292161
    KPSS p-value: 0.010000
    non-stationary - null hypothesis of stationary can be rejected at a 5% significance level

1
2

augmented_dickey_fuller_test(train['ln_rnw'])
kpss_test(train['ln_rnw'])

    ADF Statistic: -0.266563
    ADF p-value: 0.930108
    non-stationary - null hypothesis of a unit root cannot be rejected at a 5% significance level
    KPSS Statistic: 2.300497
    KPSS p-value: 0.010000
    non-stationary - null hypothesis of stationary can be rejected at a 5% significance level

None of the variables can be considered stationary with the two tests. We can therefore ask whether the variables form a cointegrated system with a given number of “common trends”. Intuitively, I would argue that there exists one cointegration equation among the three variables as in the long run, one’s amount consumed and wealth accumulated should be equal to the income that one earns.

We will discuss how to test for cointegration relationships and how to estimate them using a vector error correction model in the next blog.