<jupyter_start><jupyter_text>Forecasting Inflation Using VAR: Traditional v ML Approaches The goal of this project is to understand how feature selection from Lasso and XGBoost can improve the forecasting accuracy of vector autoregressive models (VAR). A large dataset of US economic indicators constructed by economists, <NAME> and <NAME>: https://research.stlouisfed.org/econ/mccracken/fred-databases is used. In this project, the in-sample forecasts of Consumer Price Index (CPI) is compared using three separate VARs:* VAR with variables selected by me* VAR with features selected by Lasso* VAR with features selected by XGBoost<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline 
import pandas as pd 
import seaborn as sns
sns.set()
import datetime as dt
import warnings
warnings.simplefilter('ignore')

from sklearn import (linear_model, metrics, neural_network, pipeline, preprocessing, model_selection)
import statsmodels.formula.api as sm

!pip install xgboost
from xgboost import XGBRegressor, plot_importance<jupyter_output>Requirement already satisfied: xgboost in /opt/conda/lib/python3.6/site-packages (0.82)
Requirement already satisfied: numpy in /opt/conda/lib/python3.6/site-packages (from xgboost) (1.16.2)
Requirement already satisfied: scipy in /opt/conda/lib/python3.6/site-packages (from xgboost) (1.1.0)<jupyter_text>Data visualization and exploration First, some macroeconomic indicators are visualized and explored to assist the variable selection for the VAR.<jupyter_code># Some cleaning
fred = pd.read_csv("fredmd.csv")

fred = fred.drop(0)
fred["sasdate"] = pd.to_datetime(fred["sasdate"])
fred = fred.rename(columns={"sasdate": "date"})
fred = fred.set_index("date")

# Stationarity in yearly percentage change
for c in list(fred.columns.values):
    fred[c + "_diff"] = (fred[c] - fred[c].shift(12))/fred[c].shift(12) * 100

fred.head()<jupyter_output><empty_output><jupyter_text>Checking whether the **unemployment rate**, **Fed Funds rate** and **stock market indicators** have any historical relationship with inflation.<jupyter_code>def scatter_plot(data, col_1, col_2, color):
    
    """
    Creates scatter plots for inflation, unemployment rate, Fed Funds rate, the S&P500 Composite Index
    """
    
    fig, ax = plt.subplots(figsize=(10,4))
    
    ax.scatter(x=data[col_1], y=data[col_2], alpha=0.4, color=color)
    ax.set_xlabel("Inflation", fontname="Verdana")
    
    if col_2 == "UNRATE":
        ax.set_title("Inflation and unemployment rate", fontsize=16, fontname="Verdana", loc="left")
        ax.set_ylabel("Unemployment rate", fontname="Verdana")
    elif col_2 == "S&P 500":
        ax.set_title("Inflation and S&P500", fontsize=16, fontname="Verdana", loc="left")
        ax.set_ylabel("S&P500", fontname="Verdana")
    else:
        ax.set_title("Inflation and Fed Funds rate", fontsize=16, fontname="Verdana", loc="left")
        ax.set_ylabel("Fed Funds rate", fontname="Verdana")<jupyter_output><empty_output><jupyter_text>Economists often refer to the Phillips Curve, showing a trade-off between inflation and **unemployment**. Looking at the data, is this trade-off true? Do stock markets have any predictive power for inflation? I seek to explore the relationship between the inflation rate and **S&P500** Composite index.The Fed Reserve lowers inflation by raising interest rates. Is the **Fed Funds rate** useful for prediction?<jupyter_code>for f,c in dict(zip(["UNRATE", "S&P 500", "FEDFUNDS"], ["navy", "firebrick", "orange"])).items():
    scatter_plot(data=fred, col_1="CPIAUCSL_diff", col_2=f, color=c)<jupyter_output><empty_output><jupyter_text>Despite some outliers, there appears a curved negative relation between the unemployment rate and inflation, providing some visual evidence for the Phillips Curve. This relation is broken down into various subsamples below. There's no clear relationship between stock markets and inflation, contrary to popular beliefs. This will not be too useful a variable for my VAR.There appears to be a positive upward-sloping relationship between the Fed Funds rate and inflation. This is not surprising, given that monetary policy responds to inflation (and vice versa).<jupyter_code>def scatter_subs(data, col_1, col_2, color):
    
    """
    Break down scatterplots into different years
    """
    
    data1 = data.loc["1960-01-01":"1979-12-01", :]
    data2 = data.loc["1980-01-01":"1999-12-01", :]
    data3 = data.loc["2000-01-01":"2019-01-01", :]

    fig, ax = plt.subplots(1, 3, figsize=(12,4), sharey=True)
    
    ax[0].scatter(x=data1[col_1], y=data1[col_2], alpha=0.4, color=color)
    ax[1].scatter(x=data2[col_1], y=data2[col_2], alpha=0.4, color=color)
    ax[2].scatter(x=data3[col_1], y=data3[col_2], alpha=0.4, color=color)

    ax[0].set_title("1960-1979", fontsize=14, fontname="Verdana")
    ax[1].set_title("1980-1999", fontsize=14, fontname="Verdana")
    ax[2].set_title("2000-2019", fontsize=14, fontname="Verdana")

    for i in list(range(3)):
        ax[i].set_xlabel("Inflation")
        if col_2 == "UNRATE":
            ax[i].set_ylabel("Unemployment rate (%)")
        elif col_2 == "S&P 500":
            ax[i].set_ylabel("S&P500")
        else:
            ax[i].set_ylabel("Fed funds rate")
            
scatter_subs(data=fred, col_1="CPIAUCSL_diff", col_2="UNRATE", color="navy")<jupyter_output><empty_output><jupyter_text>Outliers may be explained by trends in across 1970 and 1980. In the 70s, the US economy experienced high inflation. The clustering of low inflation values from 1980-1999 may correspond to the recession of the 1980s, following aggressive interest rate hikes by the Fed in response to the high inflation in the 70s.<jupyter_code>scatter_subs(data=fred, col_1="CPIAUCSL_diff", col_2="S&P 500", color="firebrick")<jupyter_output><empty_output><jupyter_text>Across different subsamples of years, there does not appear to be any relationship between inflation and stock markets.<jupyter_code>scatter_subs(data=fred, col_1="CPIAUCSL_diff", col_2="FEDFUNDS", color="orange")<jupyter_output><empty_output><jupyter_text>In most intervals, the Fed funds rate has a positive, upward-sloping relationship with inflation. This makes sense as interest rates rise in response to inflation. Interest rates historically have been much higher than in recent times.  VAR: Traditional approach<jupyter_code>from statsmodels.tsa.api import VAR

fred = pd.read_csv("fredmd.csv")

fred = fred.drop(0)
fred["sasdate"] = pd.to_datetime(fred["sasdate"])
fred = fred.rename(columns={"sasdate": "date"})
fred = fred.set_index("date")

# Stationarity
# First differences instead of percentage change
# It seems later on that Lasso had an issue fitting percentage change

for c in list(fred.columns.values):
    fred[c + "_diff"] = fred[c] - fred[c].shift(12)<jupyter_output><empty_output><jupyter_text>The VAR model has 6 variables: Inflation, real personal income, unemployment rate, Fed Funds rate, a trade-weighted index of the US dollar with several major currencies and total housing starts, as well as 12 lags for each variable. $B$ represents a matrix of coefficients.$$\begin{bmatrix} \pi_{t} \\ rpi_{t} \\ unemp_{t} \\ fedrate_{t} \\ usd_{t} \\ houst_{t} \end{bmatrix} = \alpha_t + B_{t-1}\begin{bmatrix} \pi_{t-1} \\ rpi_{t-1} \\ unemp_{t-1} \\ fedrate_{t-1} \\ usd_{t-1} \\ houst_{t-1} \end{bmatrix} + \dots + B_{t-12}\begin{bmatrix} \pi_{t-12} \\ rpi_{t-12} \\ unemp_{t-12} \\ fedrate_{t-12} \\ usd_{t-12} \\ houst_{t-12} \end{bmatrix}+ \epsilon_t$$<jupyter_code># Rename VAR variables

new_names = ["cpi", "rpi", "unemp", "fedrate", "usd", "houst"]
old_names = ["CPIAUCSL", "RPI", "UNRATE", "FEDFUNDS", "TWEXMMTH", "HOUST"]

for n,o in dict(zip(new_names, old_names)).items():
    fred = fred.rename(columns={o:n})
    
new_names_diff = ["cpi_diff", "rpi_diff", "unemp_diff", "fedrate_diff", "usd_diff", "houst_diff"]
old_names_diff = ["CPIAUCSL_diff", "RPI_diff", "UNRATE_diff", "FEDFUNDS_diff", "TWEXMMTH_diff", "HOUST_diff"]

for nd,od in dict(zip(new_names_diff, old_names_diff)).items():
    fred = fred.rename(columns={od:nd})
    
fred.head()
def plot_vars(data, levels, color, leveltype):
    
    """
    Displays historical trends of VAR variables
    And see if it's sensible to just select levels instead of differences
    """
    
    fig, ax = plt.subplots(1, 6, figsize=(16,2.5), sharex=True)
    
    palettes = ["blue", "green", "red", "orange", "purple", "black"]
    
    for col, i in dict(zip(levels, list(range(6)))).items():
        data[col].plot(ax=ax[i], legend=True, linewidth=1.0, color=color, sharex=True)     
    
    fig.set_facecolor("floralwhite")
    fig.suptitle(f"Historical trends of VAR {leveltype} variables", 
                 fontsize=14, fontweight="bold", fontname="Verdana")
plot_vars(fred, levels=new_names, color="royalblue", leveltype="levels")
plot_vars(fred, levels=new_names_diff, color="firebrick", leveltype="difference")<jupyter_output><empty_output><jupyter_text>Check for *stationarity*: Inflation and real personal income levels are clearly not stationary. The first differences in Fed Funds rate and trade-weighted USD index are taken. Unemployment rate in levels looks stationary.<jupyter_code>start_date = "2009-02-01"
end_date = "2019-01-01"

def var_create(columns, data):
    
    """
    Creates vector autoregressive model given data and list of selected variables
    Returns the MSE between forecasted and actual values
    Also returns the concatenated dataset for visualization purposes
    """
    
    data = data[columns]
    data = data.dropna(axis=0)
    data.index.to_period("M")
    
    # Split dataset and run VAR on the trained part
    data_train = data.loc["1973-01":"2009-01", :]
    var_train = VAR(data_train)
    results = var_train.fit(12)
    lag_order = results.k_ar
    forecasted = pd.DataFrame(results.forecast(data_train.values[-lag_order:], 120)) # Forecast 120 months
    
    # Rename forecasted columns
    forecasted_names = list(forecasted.columns.values)
    data_train_names = list(data_train.columns.values)
    
    var_dict = dict(zip(forecasted_names, data_train_names))
    
    for f,t in var_dict.items():
        forecasted = forecasted.rename(columns={f:t + "_fcast"})
        
    forecasted.index = pd.DatetimeIndex(start=start_date, periods=forecasted.shape[0], freq="MS")
    forecasted.index.names = ["date"]
    
    # Parse together forecasted data with original dataset
    final_data = pd.merge(forecasted, data, left_index=True, right_index=True)
    final_data = final_data.sort_index(axis=0, ascending=True)
    final_data = pd.concat([data_train, final_data], sort=True, axis=0)
    final_data = final_data.sort_index(axis=0, ascending=True)
    
    var_mse = metrics.mean_squared_error(final_data.loc[start_date:end_date,"cpi_diff_fcast"], 
                           final_data.loc[start_date:end_date,"cpi_diff"])
    
    return var_mse, final_data
mse1, df1 = var_create(columns=["unemp", "cpi_diff", "fedrate_diff", 
                    "usd_diff", "rpi_diff", "houst"], data=fred)

print(f"The mean squared error between the forecasted and actual values is {mse1}")
def plot_cpi(final_data, var_mse, approach):
    
    """
    Plots the actual values against forecast
    """
    
    fig, ax = plt.subplots(figsize=(14,6))
    colors = sns.color_palette("deep", 8)

    final_data["cpi_diff_fcast"].plot(ax=ax, legend=True, linewidth=2.5, linestyle="dashed")
    final_data["cpi_diff"].plot(ax=ax, legend=True, alpha=0.6, linestyle="solid")
    
    if approach=="traditional":
        ax.set_title("VAR in-sample forecast, traditional approach", fontsize=16, 
                     fontweight="bold", fontname="Verdana", loc="left")
    elif approach=="lasso":
        ax.set_title("VAR in-sample forecast, Lasso approach", fontsize=16, 
                     fontweight="bold", fontname="Verdana", loc="left")
    elif approach=="XGBoost":
        ax.set_title("VAR in-sample forecast, XGBoost approach", fontsize=16, 
                     fontweight="bold", fontname="Verdana", loc="left")
    else:
        ax.set_title("VAR in-sample forecast", fontsize=16, 
                     fontweight="bold", fontname="Verdana", loc="left")
    
    ax.set_ylabel("First differences", fontname="Verdana")
    ax.legend([f"VAR Forecast, MSE={var_mse}", "CPI Real Values"])
    
plot_cpi(final_data=df1, var_mse=mse1, approach="traditional")<jupyter_output><empty_output><jupyter_text>The MSE is high, but the visualizations show that the VAR actually performs quite well over the forecast period.  VAR: Lasso approach I now use Lasso to determine whether the features selected by VAR do a better job at predicting inflation. Does it select the same variables?<jupyter_code>def timeseries_train_test_split(X, y, testsize):
    """
    This function splits the sample into a trained and test data
    """
    
    X_train = X.loc[:"2009-01-01"] # Until a particular date
    y_train = y.loc[:"2009-01-01"]
    X_test = X.loc["2009-01-01":]
    y_test = y.loc["2009-01-01":]
    
    return X_train, y_train, X_test, y_test

def train_test_plot(model, X_train, X_test):
    """
    This will plot the actual values of CPI against the one fitted by the model
    We train the model until 2009 and then use it from 2009 onwards on the test features dataset
    """
    fig, ax = plt.subplots(figsize=(12,4))
    colors = sns.color_palette("deep", 8)
    
    yvalues = pd.DataFrame(y_test)
    
    forecasted = list(model.predict(X_test)) # Use the model fit on features data from 2009 onwards
    df_fcast = pd.DataFrame({"date": list(yvalues.index), "cpi_fcast": forecasted})
    df_fcast = df_fcast.set_index("date")
    
    df = pd.merge(yvalues, df_fcast, left_index=True, right_index=True)

    df["cpi_fcast"].plot(ax=ax, legend=True, linewidth=2.5, linestyle="dashed", color="forestgreen") # CPI fitted
    df["cpi_diff"].plot(ax=ax, legend=True, linewidth=1.5, linestyle="solid", color="salmon") # Actual CPI values
    
    ax.set_title("CPI vs. Model's CPI")
    ax.set_ylabel("First differences")
    ax.legend(["Fitted CPI","Actual CPI"])
# Split into X and y
cpi_target = fred.dropna().cpi_diff
fred_features = fred.dropna().drop(["cpi_diff"], axis=1)<jupyter_output><empty_output><jupyter_text>One thing to realize is that in time series data, Lasso performs the k-folds on a rolling basis, as the data is folded across time observations.<jupyter_code>X_train, y_train, X_test, y_test = timeseries_train_test_split(X=fred_features, y=cpi_target, testsize=0.25)

lasso = linear_model.LassoCV(cv=model_selection.TimeSeriesSplit(n_splits=5), 
                             alphas=None, tol = 10000, normalize=True) 

fred_lasso = lasso.fit(X_train, y_train)
optimal_alpha = fred_lasso.alpha_

lasso2 = linear_model.Lasso(alpha=optimal_alpha, normalize=True)
lasso2.fit(X_train, y_train)

train_test_plot(lasso2, X_train, X_test) 
metrics.mean_squared_error(y_test, lasso2.predict(X_test))<jupyter_output><empty_output><jupyter_text>The MSE is actually quite amazingly small, showing that the test data fitted into Lasso actually predicts CPI very well and closely fits with the actual values. Now for the feature selection:<jupyter_code>lasso_coefs = pd.DataFrame({"features":list(X_train), "coef": lasso2.coef_})
lasso_coefs = lasso_coefs[lasso_coefs.coef != 0.0]
lasso_coefs.sort_values("coef", ascending=False)<jupyter_output><empty_output><jupyter_text>Got the definitions from McCracken's website: https://s3.amazonaws.com/files.fred.stlouisfed.org/fred-md/Appendix_Tables_Update.pdf, but the features are renamed below. The five features that are selected by Lasso are **CPI - All items less medical**, **CPI - Services**, **CPI - Commodities**, **CPI - All items less food**, **CPI - All items less shelter** and **CPI - Medical care**. This is strange, considering concerns of multicollinearity. It's almost like a univariate regression breaking down CPI into different parts. Nonetheless, they are used in the VAR to see how they perform.<jupyter_code>old_names2 = ["CUSR0000SA0L5_diff", "CUSR0000SAS_diff", 
              "CUSR0000SAC_diff", "CPIULFSL_diff", "CUSR0000SA0L2_diff"]

new_names2 = ["cpi_lessmed", "cpi_serv", "cpi_comm", "cpi_lessfd", "cpi_lessshelt"]

for n,o in dict(zip(new_names2, old_names2)).items():
    fred = fred.rename(columns={o:n})

new_names2 = new_names2 + ["cpi_diff"]
mse2, df2 = var_create(columns=new_names2, data=fred)
print(f"The mean squared error between the forecasted and actual values is {mse2}")
plot_cpi(final_data=df2, var_mse=mse2, approach="lasso")<jupyter_output><empty_output><jupyter_text>The MSE is actually worse than that of the VAR using traditional approach. This shows that machine learning has its limits.  VAR: XGBoost approach Next, I use XGBoost to select features for the VAR. It turns out that XGBoost and other tree-based models poorly handle trends in time series datasets. Ideally, I should restrict my dataset to stationary variables so that it doesn't 'read' variables that are trending upwards.<jupyter_code>fred_features = fred_features[list(fred_features.filter(regex = "_diff"))] # Drop non-stationary variables 
X_train, y_train, X_test, y_test = timeseries_train_test_split(X=fred_features, y=cpi_target, testsize=0.25)
# Build XGBoost
# Normalize the features 
scaler = preprocessing.StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)

xgb = XGBRegressor()
xgb.fit(X_train_scaled, y_train)

X_test_scaled = scaler.transform(X_test)

train_test_plot(model=xgb, X_train=X_train_scaled, X_test=X_test_scaled)
fig_xgb, ax_xgb = plt.subplots(figsize=(10,6))
plot_importance(xgb, max_num_features=5, ax=ax_xgb)
# Feature selection

for i, name in dict(zip(list(range(len(X_train.columns))), list(X_train.columns.values))).items():
    if i == 0:
        print(f"Feature {i} is {name}")
    elif i == 99:
        print(f"Feature {i} is {name}")
    elif i == 106:
        print(f"Feature {i} is {name}")
    elif i == 109:
        print(f"Feature {i} is {name}")
    elif i == 111:
        print(f"Feature {i} is {name}")<jupyter_output>Feature 0 is rpi_diff
Feature 99 is WPSFD49207_diff
Feature 106 is CPITRNSL_diff
Feature 109 is CUSR0000SAD_diff
Feature 111 is CPIULFSL_diff<jupyter_text>It turns out that **real personal income**, **PPI - Finished Goods**, **CPI - Transportation**, and **personal consumption expenditure (PCE)**, and **PCE - Non-durable goods** are the most important features for the XGBoost Regressor. The variables are renamed for convenience.<jupyter_code>old_names3 = ["AMDMUOx", "cpi_lessmed", 
              "cpi_lessfd", "rpi", "PCEPI_diff"]

new_names3 = ["orders", "cpi_lessmed", "cpi_lessfd", "rpi", "pce_diff"]

for n,o in dict(zip(new_names3, old_names3)).items():
    fred = fred.rename(columns={o:n})

new_names3 = new_names3 + ["cpi_diff"]
mse3, df3 = var_create(columns=new_names3, data=fred)
print(f"The mean squared error between the forecasted and actual values is {mse3}")
plot_cpi(final_data=df3, var_mse=mse3, approach="XGBoost")<jupyter_output><empty_output><jupyter_text>Wow! The forecast is very poor. This highlights the challenges of using machine learning to select variables for forecasting.  Conclusions From this exercise, I learned that there's a fundamental difference between prediction/fit and forecasting. While the historical prediction of a model might be very good and training the model on a test dataset might yield low MSEs, forecasting with those same features is tricky. Because variables in VAR interact endogenously with each other, there's much uncertainty with how the selected features will trend over time. When it comes to forecasting, one needs to decide which variables are the most important and not just rely solely on fit. A summary of the MSEs between actual and forecasted values for each approach is provided below.<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ax = sns.barplot(x = ["Traditional", "Lasso", "XGB"], y = [mse1, mse2, mse3])
ax.set_title("Comparing MSEs", fontsize=16, fontname="Verdana")
ax.set_ylabel("MSE", fontname="Verdana")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12009318001586043, 'fraction_numerical': 0.02359238699444885, 'mean_word_length': 4.781375358166189, 'pattern_counts': {'":': 19, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '633827', 'license': 'MIT', 'n_tokens_mistral': 7040, 'n_tokens_neox': 6633, 'n_words': 1848, 'path': 'Vector Autoregressive Inflation Forecasting.ipynb', 'repo_name': 'atheesh1998/Time-Series-Analysis', 'size': '20190'}