<jupyter_start><jupyter_text>______Copyright <NAME>For more information, visit us at www.pieriandata.com  VARMA(p,q) Vector Autoregressive Moving AverageThis lesson picks up where VAR(p) left off.Recall that the system of equations for a 2-dimensional VAR(1) model is:&nbsp;&nbsp;&nbsp;&nbsp;$y_{1,t} = c_1 + \phi_{11,1}y_{1,t-1} + \phi_{12,1}y_{2,t-1} + \varepsilon_{1,t}$&nbsp;&nbsp;&nbsp;&nbsp;$y_{2,t} = c_2 + \phi_{21,1}y_{1,t-1} + \phi_{22,1}y_{2,t-1} + \varepsilon_{2,t}$where the coefficient $\phi_{ii,l}$ captures the influence of the $l$th lag of variable $y_i$ on itself,the coefficient $\phi_{ij,l}$ captures the influence of the $l$th lag of variable $y_j$ on $y_i$.Most importantly, $\varepsilon_{1,t}$ and $\varepsilon_{2,t}$ are white noise processes that may be correlated.In a VARMA(p,q) model we give the error terms $\varepsilon_t$ a moving average representation of order $q$. FormulationWe've seen that an autoregressive moving average ARMA(p,q) model is described by the following:&nbsp;&nbsp;&nbsp;&nbsp; $y_{t} = c + \phi_{1}y_{t-1} + \phi_{2}y_{t-2} + \dots + \phi_{p}y_{t-p} + \theta_{1}\varepsilon_{t-1} + \theta_{2}\varepsilon_{t-2} + \dots + \theta_{q}\varepsilon_{t-q} + \varepsilon_{t}$ A $K$-dimensional VARMA model of order $(p,q)$ considers each variable $y_K$ in the system.For example, the system of equations for a 2-dimensional VARMA(1,1) model is:&nbsp;&nbsp;&nbsp;&nbsp;$y_{1,t} = c_1 + \phi_{11,1}y_{1,t-1} + \phi_{12,1}y_{2,t-1} + \theta_{11,1}\varepsilon_{1,t-1} + \theta_{12,1}\varepsilon_{2,t-1} + \varepsilon_{1,t}$&nbsp;&nbsp;&nbsp;&nbsp;$y_{2,t} = c_2 + \phi_{21,1}y_{1,t-1} + \phi_{22,1}y_{2,t-1} + \theta_{21,1}\varepsilon_{1,t-1} + \theta_{22,1}\varepsilon_{2,t-1} + \varepsilon_{2,t}$where the coefficient $\theta_{ii,l}$ captures the influence of the $l$th lag of error $\varepsilon_i$ on itself,the coefficient $\theta_{ij,l}$ captures the influence of the $l$th lag of error $\varepsilon_j$ on $\varepsilon_i$,and $\varepsilon_{1,t}$ and $\varepsilon_{2,t}$ are residual white noise.The general steps involved in building a VARMA model are:* Examine the data* Visualize the data* Test for stationarity* If necessary, transform the data to make it stationary* Select the appropriate p and q orders* Instantiate the model and fit it to a training set* If necessary, invert the earlier transformation* Evaluate model predictions against a known test set* Forecast the future Related Functions:varmax.VARMAX(endog[, exog, order, trend, …])&nbsp;&nbsp;Vector Autoregressive Moving Average with eXogenous regressors modelvarmax.VARMAXResults(model, params[, …])&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Class to hold results from fitting an VARMAX model   For Further Reading:Statsmodels Tutorial:&nbsp;&nbsp;Time Series Analysis by State Space MethodsStatsmodels Example:&nbsp;&nbsp;VARMAX models  Perform standard imports and load datasetFor this analysis we'll reuse our money and spending datasets. We'll look at the M2 Money Stock which is a measure of U.S. personal assets, and U.S. personal spending. Both datasets are in billions of dollars, monthly, seasonally adjusted. They span the 21 years from January 1995 to December 2015 (252 records).Sources: https://fred.stlouisfed.org/series/M2SL https://fred.stlouisfed.org/series/PCE<jupyter_code>import numpy as np
import pandas as pd
%matplotlib inline

# Load specific forecasting tools
from statsmodels.tsa.statespace.varmax import VARMAX, VARMAXResults
from statsmodels.tsa.stattools import adfuller
from pmdarima import auto_arima
from statsmodels.tools.eval_measures import rmse

# Ignore harmless warnings
import warnings
warnings.filterwarnings("ignore")

# Load datasets
df = pd.read_csv('../UnusedData/M2SLMoneyStock.csv',index_col=0, parse_dates=True)
df.index.freq = 'MS'

sp = pd.read_csv('../UnusedData/PCEPersonalSpending.csv',index_col=0, parse_dates=True)
sp.index.freq = 'MS'<jupyter_output><empty_output><jupyter_text>Inspect the data<jupyter_code>df = df.join(sp)
df.head()
df = df.dropna()
df.shape<jupyter_output><empty_output><jupyter_text>Plot the source data<jupyter_code>title = 'M2 Money Stock vs. Personal Consumption Expenditures'
ylabel='Billions of dollars'
xlabel=''

ax = df['Spending'].plot(figsize=(12,5),title=title,legend=True)
ax.autoscale(axis='x',tight=True)
ax.set(xlabel=xlabel, ylabel=ylabel)
df['Money'].plot(legend=True);<jupyter_output><empty_output><jupyter_text>Test for stationarity, perform any necessary transformationsIn the previous section we applied the augmented Dickey-Fuller test and found that a second-order difference achieved stationarity. In this section we'll perform the auto_arima prediction to identify optimal $p$ and $q$ orders.<jupyter_code># INCLUDED HERE IF YOU CHOOSE TO USE IT
def adf_test(series,title=''):
    """
    Pass in a time series and an optional title, returns an ADF report
    """
    print(f'Augmented Dickey-Fuller Test: {title}')
    result = adfuller(series.dropna(),autolag='AIC') # .dropna() handles differenced data
    
    labels = ['ADF test statistic','p-value','# lags used','# observations']
    out = pd.Series(result[0:4],index=labels)

    for key,val in result[4].items():
        out[f'critical value ({key})']=val
        
    print(out.to_string())          # .to_string() removes the line "dtype: float64"
    
    if result[1] <= 0.05:
        print("Strong evidence against the null hypothesis")
        print("Reject the null hypothesis")
        print("Data has no unit root and is stationary")
    else:
        print("Weak evidence against the null hypothesis")
        print("Fail to reject the null hypothesis")
        print("Data has a unit root and is non-stationary")<jupyter_output><empty_output><jupyter_text>NOTE:  When performing the auto_arima function we're likely to see a ConvergenceWarning: Maximum Likelihood optimization failed to converge. This is not unusual in models which have to estimate a large number of parameters. However, we can override the maximum iterations default of 50, and pass an arbitrarily large number with maxiter=1000. We'll see this come up again when we fit our model.<jupyter_code>auto_arima(df['Money'],maxiter=1000)
auto_arima(df['Spending'],maxiter=1000)<jupyter_output><empty_output><jupyter_text>It looks like a VARMA(1,2) model is recommended. Note that the $d$ term (2 for Money, 1 for Spending) is about to be addressed by transforming the data to make it stationary. As before we'll apply a second order difference.<jupyter_code>df_transformed = df.diff().diff()
df_transformed = df_transformed.dropna()
df_transformed.head()
len(df_transformed)<jupyter_output><empty_output><jupyter_text>Train/test/splitIt is useful to define a number of observations variable for our test set. For this analysis, let's use 12 months.<jupyter_code>nobs=12
train, test = df_transformed[0:-nobs], df_transformed[-nobs:]
print(train.shape)
print(test.shape)<jupyter_output>(238, 2)
(12, 2)<jupyter_text>Fit the VARMA(1,2) ModelThis may take awhile.<jupyter_code>model = VARMAX(train, order=(1,2), trend='c')
results = model.fit(maxiter=1000, disp=False)
results.summary()<jupyter_output><empty_output><jupyter_text>Predict the next 12 valuesUnlike the VAR model we used in the previous section, the VARMAX .forecast() function won't require that we pass in a number of previous observations, and it will provide an extended DateTime index.<jupyter_code>df_forecast = results.forecast(12)
df_forecast<jupyter_output><empty_output><jupyter_text>Invert the TransformationRemember that the forecasted values represent second-order differences. To compare them to the original data we have to roll back each difference. To roll back a first-order difference we take the most recent value on the training side of the original series, and add it to a cumulative sum of forecasted values. When working with second-order differences we first must perform this operation on the most recent first-order difference.Here we'll use the nobs variable we defined during the train/test/split step.<jupyter_code># Add the most recent first difference from the training side of the original dataset to the forecast cumulative sum
df_forecast['Money1d'] = (df['Money'].iloc[-nobs-1]-df['Money'].iloc[-nobs-2]) + df_forecast['Money'].cumsum()

# Now build the forecast values from the first difference set
df_forecast['MoneyForecast'] = df['Money'].iloc[-nobs-1] + df_forecast['Money'].cumsum()
# Add the most recent first difference from the training side of the original dataset to the forecast cumulative sum
df_forecast['Spending1d'] = (df['Spending'].iloc[-nobs-1]-df['Spending'].iloc[-nobs-2]) + df_forecast['Spending'].cumsum()

# Now build the forecast values from the first difference set
df_forecast['SpendingForecast'] = df['Spending'].iloc[-nobs-1] + df_forecast['Spending'].cumsum()
df_forecast
pd.concat([df.iloc[-12:],df_forecast[['MoneyForecast','SpendingForecast']]],axis=1)<jupyter_output><empty_output><jupyter_text>Plot the results<jupyter_code>df['Money'][-nobs:].plot(figsize=(12,5),legend=True).autoscale(axis='x',tight=True)
df_forecast['MoneyForecast'].plot(legend=True);
df['Spending'][-nobs:].plot(figsize=(12,5),legend=True).autoscale(axis='x',tight=True)
df_forecast['SpendingForecast'].plot(legend=True);<jupyter_output><empty_output><jupyter_text>Evaluate the model&nbsp;&nbsp;&nbsp;&nbsp;$RMSE = \sqrt{{\frac 1 L} \sum\limits_{l=1}^L (y_{T+l} - \hat y_{T+l})^2}$where $T$ is the last observation period and $l$ is the lag.<jupyter_code>RMSE1 = rmse(df['Money'][-nobs:], df_forecast['MoneyForecast'])
print(f'Money VAR(5) RMSE: {RMSE1:.3f}')
RMSE2 = rmse(df['Spending'][-nobs:], df_forecast['SpendingForecast'])
print(f'Spending VAR(5) RMSE: {RMSE2:.3f}')<jupyter_output>Spending VAR(5) RMSE: 243.777<jupyter_text>Clearly these results are less accurate than our earlier VAR(5) model. Still, this tells us something! Let's compare these results to individual ARMA(1,2) models<jupyter_code>from statsmodels.tsa.arima_model import ARMA,ARMAResults<jupyter_output><empty_output><jupyter_text>Money<jupyter_code>model = ARMA(train['Money'],order=(1,2))
results = model.fit()
results.summary()
start=len(train)
end=len(train)+len(test)-1
z1 = results.predict(start=start, end=end).rename('Money')
z1 = pd.DataFrame(z1)
z1<jupyter_output><empty_output><jupyter_text>Invert the Transformation, Evaluate the Forecast<jupyter_code># Add the most recent first difference from the training set to the forecast cumulative sum
z1['Money1d'] = (df['Money'].iloc[-nobs-1]-df['Money'].iloc[-nobs-2]) + z1['Money'].cumsum()

# Now build the forecast values from the first difference set
z1['MoneyForecast'] = df['Money'].iloc[-nobs-1] + z1['Money1d'].cumsum()
z1
RMSE3 = rmse(df['Money'][-nobs:], z1['MoneyForecast'])

print(f'Money VARMA(1,2) RMSE: {RMSE1:.3f}')
print(f'Money  ARMA(1,2) RMSE: {RMSE3:.3f}')<jupyter_output>Money VARMA(1,2) RMSE: 422.942
Money  ARMA(1,2) RMSE: 32.236<jupyter_text>Personal Spending<jupyter_code>model = ARMA(train['Spending'],order=(1,2))
results = model.fit()
results.summary()
start=len(train)
end=len(train)+len(test)-1
z2 = results.predict(start=start, end=end).rename('Spending')
z2 = pd.DataFrame(z2)
z2<jupyter_output><empty_output><jupyter_text>Invert the Transformation, Evaluate the Forecast<jupyter_code># Add the most recent first difference from the training set to the forecast cumulative sum
z2['Spending1d'] = (df['Spending'].iloc[-nobs-1]-df['Spending'].iloc[-nobs-2]) + z2['Spending'].cumsum()

# Now build the forecast values from the first difference set
z2['SpendingForecast'] = df['Spending'].iloc[-nobs-1] + z2['Spending1d'].cumsum()
z2
RMSE4 = rmse(df['Spending'][-nobs:], z2['SpendingForecast'])

print(f'Spending VARMA(1,2) RMSE: {RMSE2:.3f}')
print(f'Spending  ARMA(1,2) RMSE: {RMSE4:.3f}')<jupyter_output>Spending VARMA(1,2) RMSE: 243.777
Spending  ARMA(1,2) RMSE: 52.334
{'fraction_non_alphanumeric': 0.13523520996381386, 'fraction_numerical': 0.023899688630817133, 'mean_word_length': 6.363073110285006, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 66, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, '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': '269692', 'license': 'MIT', 'n_tokens_mistral': 4159, 'n_tokens_neox': 3912, 'n_words': 1101, 'path': 'tsa/jose/TSA_COURSE_NOTEBOOKS/06-General-Forecasting-Models/09-Vector-AutoRegressive-Moving-Average-VARMA.ipynb', 'repo_name': 'juspreet51/ml_templates', 'size': '11889'}