Forecasting retail turnover of household goods in the Australian economy
Data taken from the ABS website is forecasted 4 periods into the future by different models
- Ian Dias
- Nov 2020
- Introduction
- 1. Visualization of the time-series
- Visualizing the data and patterns
- 2. Stationarity checks
- 3. Predictors
- 3.1 - The benchmark: Naive Predictor
- 3.2 - AR(5, 12) Predictor
- 3.3 Linear Regression - Seasonality Controls
- 3.4 Holt Winters Smoothing
- 3.5 Adding external variables
- 3.6 Vector AutoRegression
- 4. Summary
- 5. Predictions
- Point forecasts
- Density Forecasts
- One-step ahead forecast for evaluation
I plan to forecast the monthly retail turnover of household goods in the Australian wide economy. The ABS has collected this data from May of 1982 and the last data available is until August 2020. The amount of household goods being traded is an important indication of an economy - as the goods range from electrical, to househod building and furniture. It appears that the more individuals are able to spend on these goods the seemingly more satisfaction and utility they achieve.
Due to the recent impacts of the COVID-19 pandemic, it would be of interest how this has affected the household goods industry. The Australian economy as a whole has reduced, but given that consumers are spending more time at home the impact it has had on this specific industry might not be completely obvious.
The data is accessed from: https://www.abs.gov.au/statistics/industry/retail-and-wholesale-trade/retail-trade-australia/latest-release
import pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.graphics.tsaplots import plot_pacf
sns.set_style('darkgrid')
plt.rcParams["figure.figsize"] = (10,5)
def reset(data='850101.xls'):
df = pd.read_excel(data,usecols=[0,2],sheet_name = 'Data1',skiprows='5')
df = df[10:]
df.set_index('Unnamed: 0', inplace=True)
df.index.name = 'Dates'
df.columns = ['Household goods']
df.index = pd.to_datetime(df.index)
return df
df = reset()
def plot(x,title=None,baseline=None, **kwargs):
plt.plot(x, **kwargs)
if baseline is not None:
plt.plot(baseline, label='Baseline', alpha=0.5)
plt.legend()
plt.title(title)
plot(df, title='Retail Household Goods')
There seems to be both clear seasonality and and upwards trend to the data. The data will be isolated by year to identify the seasonality patterns and the impacts of the trend.
df['year'] = df.index.year
df['month'] = df.index.month
years = df['year'].unique()
for i, y in enumerate(years):
plt.plot('month', 'Household goods', data=df[df['year'] == y], label=y)
if i % 4 == 0 or y == 2020:
plt.text(df[df['year'] == y].shape[0], df[df['year'] == y]['Household goods'][-1:].values[0], y)
Looking at the above graph, it appears that there is a definite upward tick for December and there could possibly be monthly or quarterly variation in the data.The mean also appears to be higher every year. There seems to be a significant increase in spending in 2020 starting from February and going on until the last value of October. This is most likely due to the COVID-19 virus that is impacting the world and began its Australian penetration in March. As more individuals are under lockdown they are staying home longer and this may increase the quantity of household goods they would typically require.
df = reset()
from statsmodels.tsa.stattools import adfuller
first_diff = df.diff()
def dickey_fuller(x):
'''Dickey-fuller test to check for stationarity'''
result = adfuller(x)
print('ADF Test statistic: %f' % result[0])
print('p-value: %f' % result[1])
print('The critical values are: {}' .format(result[4]))
for key,value in result[4].items():
if value > result[0]:
print('There is a unit root present and the null hypothesis is rejected. The series is stationary.')
return
print('The null hypothesis is failed to be rejected. The series is not stationary.')
dickey_fuller(df)
As the series is not stationary, the first difference will be taken to remove the impact of the trend.
dickey_fuller(df.diff().dropna())
Even after the first difference, the series is not stationarity. The log will first be taken, and then the first difference.
np_df = np.array(df.dropna(), dtype='float').squeeze()
log_df = np.log(np_df)
log_diff_df = np.diff(log_df)
dickey_fuller(log_diff_df)
This gives a unit root and the series is stationary. Taking the log helps stablise the variance and differences help stabilize the mean.
plt.plot(log_diff_df)
plt.title('First Difference of the Log of Household Goods')
Correlation functions will be used to determine the appropriate lag length.
def acf(x):
'''Plotting the correlation of lag k with lag 0'''
N = x.shape[0]
ac = x-x.mean()
lag0 = np.dot(ac, ac)
lagk = np.correlate(ac, ac, mode='full')[N-1:]
r = lagk/lag0
return r
plt.title('ACF')
lags = 20
plt.scatter(x = range(0,lags+1), y=acf(log_diff_df)[:lags+1])
Using acf given by statsmodels we get the following prettier graph:
fig, axs = plt.subplots()
plot_acf(log_diff_df, ax=axs, lags=20)
plt.show()
From the ACF above, lags 5, 12 have significant correlation with the current value and as well will be included in any further lagged process. There doesn't seem to be a geometric decay (only for past values of itself) and so it make indicate an MA process.
plot_pacf(log_diff_df, lags=20)
plt.show()
If the series was an MA process, there would be a geometric decay in the PACF.
df = reset()
df = df[1:]
df['log'] = log_diff_df
df = df.drop('Household goods',axis=1)
df = df.rename(columns={'log':'Household goods'})
The naive predictor used will have the following model specification: $y_t$ = $y_{t-1}$ and will be the baseline for any future model comparisons
plot(df, label = 'Naive',baseline=df.shift(1), title='Retail Household Goods')
def MSFE_base(yhat, model, h=1, df=df, printf=True):
MSFE = np.mean(np.square(yhat - df['Household goods'][h:]))
if printf:
return 'The MSFE of this {} is {:2f}'.format(model, MSFE)
else:
return MSFE
MSFE_base(df.shift(1)['Household goods'], 'Naive Predictor')
msfe_naive = MSFE_base(df.shift(1)['Household goods'], model = 'naive',printf=False)
As established in the ACF and PACF above, lags 5 and 12 will be used within an AR model.
$ y_t = \beta_1y_{t-5} + \beta_2y_{t-12} + \epsilon_t$
df = reset()
df = df[1:]
df['log'] = log_diff_df
df = df.drop('Household goods',axis=1)
df = df.rename(columns={'log':'Household goods'})
ytminus5 = df.shift(5)
ytminus12 = df.shift(12)
ytminus24 = df.shift(24)
df = df.join(ytminus5, rsuffix='yt-5')
df = df.join(ytminus12, rsuffix='yt-12')
df = df.dropna()
y = df['Household goods']
X = df.drop('Household goods', axis=1)
def MSFE(X, y, T0, h, plot=False, model=None):
T = len(y)
syhat = []
for i in range(T0, T-h):
temp_y = y[:i]
temp_X = X[:i]
temp_X = np.array(temp_X, dtype='float')
beta_hat = (np.linalg.inv((temp_X.T@temp_X)) @ temp_X.T@temp_y)
y_hat2 = temp_X @ beta_hat
syhat.append(y_hat2[-1])
if plot:
return syhat
return np.mean(np.square(y[T0+h:]-syhat))
def dummies(quarter = False, month = False):
X = df.reset_index()
if quarter:
X['Dates'] = X['Dates'].dt.quarter
if month:
X['Dates'] = X['Dates'].dt.month
y = X['Household goods']
X = X.drop('Household goods', axis=1)
X = pd.get_dummies(X['Dates'])
X.reset_index(inplace=True)
return X
def testing(X, T0=50, h=1, df = df):
syhat = []
for t in range(T0, len(df)-h):
y = np.array(df['Household goods'][1:t])
xt = X[1:t]
beta_hat = np.linalg.inv((xt.T@xt)) @ xt.T@y
#Now need a y_hat of t+h
y_hat = X.iloc[t+h]@np.array(beta_hat)
syhat.append(y_hat)
return syhat
def adding_index(df, yhat, T0=50, h=1):
df = df.reset_index()
df = df.join(pd.DataFrame(yhat, index=range(T0+h,len(df))))
df.set_index(df['Dates'], inplace=True)
df = df.rename(columns={0:'yhat'})
return df
def MSFE_new(syhat, model, T0=50, h=1):
syhat = np.array(syhat)
ytph = df['Household goods'][T0+h]
msfe = np.mean(np.square(ytph-syhat))
return 'The MSFE of {} is {}'.format(model, msfe)
AR = MSFE(X, y, 50, 1, plot=True)
df_new = adding_index(df, AR)
df_new = df_new.dropna()
plot(df_new['yhat'], baseline=df_new['Household goods'], label='AR Model')
model = 'Auto Regressive Model'
print('The MSFE of this {} is {:2f}'.format(model, MSFE(X,y,50,1)))
msfe_ar = MSFE(X,y,50,1)
Two model specifications will be used, one to act as a dummy for yearly quarters and the other for the months.
$ S_1: y_t = a_1 + a_1t + \sum_{i=1}^{4}\alpha_iD_{it} + \epsilon_t $
$ S_2: y_t = a_1t + \sum_{i=1}^{12}\beta_iM_{it}+\epsilon_t $
$S_1$: This includes a trend - $t$, and a dummy for all quarters $Di$, where $i$ $\epsilon$ [1,4]
$S_2$: This includes a trend $t$, and a dummy for all months of the year - $M_i$, where $i$ $\epsilon$ [1,12]
X_2 = dummies(quarter=True, month=False)
syhat_s2 = testing(X_2)
plot(adding_index(df, syhat_s2)['yhat'], baseline=df['Household goods'], label='First Specification')
MSFE_base(syhat_s2, 'First Specification',df=df[50:])
msfe_s1 = MSFE_base(syhat_s2, 'frs', df = df[50:],printf=False)
X_4 = dummies(quarter=False, month=True)
syhat_s4 = testing(X_4)
plot(adding_index(df, syhat_s4)['yhat'], baseline=df['Household goods'], label='Second Specification')
MSFE_base(syhat_s4, 'Second Specification',df=df[50:])
msfe_s2 = MSFE_base(syhat_s4, 'frs', df = df[50:],printf=False)
def HWS(T0, s, alpha, beta, gamma, h, T = len(y)):
syhat = []
S = np.zeros((1,T-h))[0]
L = np.mean(y[1:s])
b = 0
S[1:s] = y[1:s] - L
for t in range(s+1, T-h):
newL = alpha * (y[t] - S[t-s]) + (1-alpha) * (L+b)
newb = beta * (newL-L) + (1-beta)*b
S[t] = gamma * (y[t] - newL) + (1-gamma)*S[t-s]
yhat = newL + h*newb + S[t+h-s]
L = newL
b = newb
if t >= T0:
syhat.append(yhat)
return syhat
y = df['Household goods']
hws_syhat = HWS(50, 4, 0.4, 0.4, 0.4, 1)
#df_new = adding_index(df, hws_syhat, T0=64)
df_new = adding_index(df, hws_syhat, T0=50)
df_new = df_new.dropna()
df_new = df_new.drop('Dates', axis=1)
plot(df_new['yhat'], baseline=df['Household goods'], label='Holt Winters Smoothing')
MSFE_base(df_new['yhat'], 'Holt Winters Smoothing',df=df)
def optim_point(T0):
alpha = np.arange(0,1,0.1)
h = 1
beta = gamma = alpha
all_values = []
index = []
temp = []
for i in alpha:
for j in beta:
for k in gamma:
syhat = HWS(50, 4, i, j, k, h)
temp = np.mean(np.square(syhat - df['Household goods'][T0+h:]))
index.append((i,j,k))
all_values.append(temp)
opt_point = np.array(index)[all_values == np.min(all_values)]
return opt_point
optimal = optim_point(50)
print('The optimal parameters for the Holt Winter Smoothing are: \n{}'.format(optimal))
hws_syhat = HWS(T0=50, s=4, alpha=0, beta=0, gamma=0.1, h=1)
df_new = adding_index(df, hws_syhat, T0=50)
df_new = df_new.dropna()
plot(df_new['yhat'], baseline=df['Household goods'], label='Holt Winters Smoothing')
msfe_hw = MSFE_base(df_new['yhat'], 'Holt Winters Smoothing',df=df, printf=False)
MSFE_base(df_new['yhat'], 'Holt Winters Optimized Smoothing',df=df)
The Australian National Accounts: National Income, Expenditure and Product dataset was downloaded from the ABS through the following link: https://www.abs.gov.au/statistics/economy/national-accounts/australian-national-accounts-national-income-expenditure-and-product/latest-release#data-download
It gives information of a broader economic glow in Australia. As this heavily relate to the amount of Household Goods being bought, it was thought that adding these as regressors would aid any future predictions of Household Goods. However, there are 108 potential regressors within this new dataset. Principal Component Analysis (PCA) along with relevant Economic knowledge will be used to determine which regressors to add
df_nat = pd.read_excel('National_accounts.xls',sheet_name = 'Data1', index_col=0, prase_dates=True)
Xs = df_nat[10:].dropna(axis=1)
yXs = Xs.join(df).dropna()
new = (yXs).drop('Household goods', axis=1)
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
sc = StandardScaler()
X_std = sc.fit_transform(new)
pca = PCA(n_components=0.95)
X_pca = pca.fit_transform(X_std)
n_pcs= pca.n_components_ # get number of component
# get the index of the most important feature on EACH component
most_important = [np.abs(pca.components_[i]).argmax() for i in range(n_pcs)]
initial_feature_names = new.columns
# get the most important feature names
most_important_names = [initial_feature_names[most_important[i]] for i in range(n_pcs)]
Through the above PCA analysis, the most relevant columns are the ones below.
most_important_names
The new timeseries have values collected only every 3 months while household goods have values collected every month. Instead of reducing the number of household good observations to match the new timeseries, it was decided to instead interpolate the new timeseries datasets to fill in the gaps. The method used was a simple linear method as it was determined the underlying functions of RGDI, Terms of Trade and Household saving ratio were not too ambitious.
reg = yXs[['RGDI', 'Terms of trade: Index - Percentage changes ;.1', 'Household saving ratio: Ratio ;.2']]
y = df['Household goods']
all_df = df.join(reg)
all_df = all_df.apply(pd.to_numeric)
all_df = all_df.interpolate(method='linear', limit_direction='forward', axis=0)[1:]
all_df = all_df.rename(columns={'Terms of trade: Index - Percentage changes ;.1':'Terms of trade',
'Household saving ratio: Ratio ;.2':'Savings ratio'})
new_X = all_df[['RGDI','Terms of trade','Savings ratio']]
new_y = pd.DataFrame(all_df['Household goods'])
Equation 1: $ y_t = \beta_0 + \beta_1RGDI + \beta_2TOT + \beta_3Savings $
var_syhat = testing(new_X, df=new_y)
plot(adding_index(all_df, var_syhat)['yhat'], title='Equation 1',baseline=(all_df['Household goods'][51:]))
Equation 2: $ y_t = y_{t-1} \beta_0 + \beta_1RGDI + \beta_2TOT + \beta_3Savings $
AR_X = new_X.join(y.shift(1))
AR_syhat = testing(AR_X, df=new_y)
plot(adding_index(all_df, AR_syhat)['yhat'], title='Equation 2', baseline=(all_df['Household goods']))
msfe_ar1 = MSFE_base(AR_syhat, 'Linear Regression',df=df[51:], printf=False)
MSFE_base(AR_syhat, 'Linear Regression',df=df[51:])
Equation 3: $ y_t = \alpha_0 y_{t-1} + \alpha_1RGDI + \alpha_2TOT + \alpha_3Savings + a_1t + \sum_{i=1}^{12} \beta_iM_{it} + \epsilon_t$
X_4 = X_4[1:].set_index(all_df.index)
new_X_4 = all_df.join(X_4)
new_X_4[['Household goods', 'Household goodsyt-5', 'Household goodsyt-12']] = new_X_4[['Household goods', 'Household goodsyt-5', 'Household goodsyt-12']].shift(1)
new_X_4 = new_X_4[1:]
new_X_4_syhat = testing(new_X_4, df=new_y[1:])
len(new_X_4_syhat)
plot(adding_index(all_df[1:], new_X_4_syhat)['yhat'], title = 'Equation 3', baseline=(all_df['Household goods'][1:]))
msfe_ar2 = MSFE_base(new_X_4_syhat, 'Equation 3',df=df[52:],printf=False)
MSFE_base(new_X_4_syhat, 'Equation 3',df=df[52:])
dickey_fuller(all_df['RGDI'])
dickey_fuller(all_df['Terms of trade'])
dickey_fuller(all_df['Savings ratio'])
All relevant variables [RGDI, Terms of Trade, Savings ratio] are stationary.
$ y_t = b + B_1y_{t-1} + ... + B_py_{t-p} + \epsilon_t$
where $b$ is an n x 1 vector of intercepts and $B_i$ is an n x n coefficient matrix
Rewrite our model as $ y = X\beta + \epsilon $, where $\epsilon ~ N(0,\Sigma)$
$ \beta = (X'X)^{-1}X'y $
$ \Sigma = 1/T(y_t - X_t\beta)(y_t - X_t\beta)' $
X = all_df[['Household goods', 'RGDI', 'Terms of trade', 'Savings ratio']]
X['Household goods'] = X['Household goods'].shift(1)
X = X.rename(columns={'Household goods': 'Household goods -1'})
X = X[1:]
y = all_df['Household goods'][1:]
Betahat = np.linalg.inv(X.T@X) @ X.T@y
Betahat = np.array(Betahat).reshape(-1,1)
T = len(y)
T0 = 50
yhatVAR = []
for t in range(T0, T):
yt = y[:t]
Xt = X[:t]
beta = np.array(np.linalg.inv(X.T@X) @ X.T@y).reshape(-1,4)
yt = np.array(yt)
yt = yt.reshape(-1,1)
all_yt = np.concatenate((yt[:t-3], yt[1:t-2], yt[2:t-1], yt[3:t]), axis=1, ).reshape(-1,4)
yhatVAR.append((beta @ all_yt.T)[-1][-1])
msfe_var = np.mean(np.square(y[T0+1:] - yhatVAR[:-1]))
print('The MSFE of the VAR is: ', msfe_var)
y = pd.DataFrame(y)
y.loc['2020-09-01 00:00:00'] = yhatVAR[-1]
yhatVARindex = adding_index(df=y, yhat = yhatVAR)['yhat']
plot(yhatVARindex, baseline=df['Household goods'], title='Vector Autoregression')
print('''Summary of MSFEs from all the models:\n\nThe Naive model: {:3f}\nThe AR(5,12) model: {:3f}\nThe Holt winters model: {:3f}
The AR and external variable (1) model: {:3f}\nThe AR and extgenal variable (2) model: {:3f}
The VAR model: {:3f}'''.format(msfe_naive, msfe_ar, msfe_hw, msfe_ar1, msfe_ar2, msfe_var))
The model with the lowest MSFE was the following:
$ y_t = \alpha_0 y_{t-1} + \alpha_1RGDI + \alpha_2TOT + \alpha_3Savings + a_1t + \sum_{i=1}^{12} \beta_iM_{it} + \epsilon_t$
def testing(X, T0=50, h=1, df = df):
syhat = []
for t in range(T0, len(df)-h):
y = np.array(df['Household goods'][1:t])
xt = X[1:t]
beta_hat = np.linalg.inv((xt.T@xt)) @ xt.T@y
#Now need a y_hat of t+h
y_hat = X.iloc[t+h]@np.array(beta_hat)
syhat.append(y_hat)
return syhat
new_X_4_syhat = testing(new_X_4, df=new_y[1:], h=1)
msfe_ar2 = MSFE_base(new_X_4_syhat, 'model',df=df[52:],printf=False)
MSFE_base(new_X_4_syhat, 'model',df=df[52:])
date_index = ['2020-09-01 00:00:00', '2020-10-01 00:00:00', '2020-11-01 00:00:00', '2020-12-01 00:00:00']
h = 4
new_X_4_syhat = testing(new_X_4, df=new_y[1:], h=h)
new_x = adding_index(all_df[8:], new_X_4_syhat[:-h])
MSFE_base(new_X_4_syhat[:-h], 'model for h=4',df=df[59:]) #MSFE for when h = 4
h=4
forecast = pd.DataFrame(new_X_4_syhat[-h:],index=date_index[:h])
forecast = forecast.rename(columns={0:'yhat'})
forecast = new_x.append(forecast)
forecast.index = pd.to_datetime(forecast.index)
plot(forecast['yhat'], title='Forecast for remaining 2020', baseline=df['Household goods'], label='forecast')
Forecasting 4 periods in the future for the remaining months of 2020 yields the following:
forecast[-h:]['yhat']
Now we need to convert it back to the original data values
orig = reset() #This is the original dataset i.e. df['Household goods']
#Need to make sure the indexes match the predicted yhat, had to use orig[18:]
#Forecast['yhat'] is the predicted yhat. Make sure the indexes are correct
real_forecast = np.exp(forecast['yhat']).cumsum() + orig[18:].set_index(forecast.index)['Household goods']
plot(real_forecast, title = 'Predicted forecast', label='forecast',baseline=orig)
The forecast with the original values for the remaining months of 2020
real_forecast[-h:]
k = 6
orig = reset()
var = 1/(len(orig)-k) * sum(np.square((real_forecast.dropna() - orig[69:]['Household goods']).dropna()))
real_density = []
for i in real_forecast[-h:]:
real_density.append((i + 1.96*np.sqrt(var), i - 1.96*np.sqrt(var)))
print('''On 2020-09-01, there is a 95% probability of the data lying between {} and {}
On 2020-10-01, there is a 95% probability of the data lying between {} and {}
On 2020-11-01, there is a 95% probability of the data lying between {} and {}
On 2020-12-01, there is a 95% probability of the data lying between {} and {}'''
.format(real_density[0][1],real_density[0][0], real_density[1][1],real_density[1][0],
real_density[2][1],real_density[2][0], real_density[3][1],real_density[3][0]))
The following is a one-step ahead forecast, instead of a 4-step forecast for the entire year.
orig = reset()
h=1
forecast = pd.DataFrame(new_X_4_syhat[-h:],index=date_index[:h])
forecast = forecast.rename(columns={0:'yhat'})
forecast = new_x.append(forecast)
forecast.index = pd.to_datetime(forecast.index)
orig = reset()
real_forecast = np.exp(forecast['yhat']).cumsum() + orig[21:].set_index(forecast.index)['Household goods']
real_forecast[-h:]
k = 6
orig = reset()
var = 1/(len(orig)-k) * sum(np.square((real_forecast.dropna() - orig[69:]['Household goods']).dropna()))
real_density = []
for i in real_forecast[-h:]:
real_density.append((i + 1.96*np.sqrt(var), i - 1.96*np.sqrt(var)))
print('On 2020-09-01, there is a 95% probability of the data lying between {} and {}'
.format(real_density[0][1],real_density[0][0]))
reset('new_data.xls')