Holt-Winters analysis
Abstract
In this post, I am going to analyse the underlying structure of the Nottingham castle temperature collected between 1920 and 1939 already analysed in previous post. In the specific my intention is to build a Holt-Winters model of the original time series. I will then analyse the remainder time series to spot further insights in its structure.
Analysis
First, I load the corresponding dataset and create a time series object from. Also defining the function to be used for running some basic statistics tests.
suppressPackageStartupMessages(library(FitARMA))
suppressPackageStartupMessages(library(forecast))
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(timeSeries))
suppressPackageStartupMessages(library(TSA))
suppressPackageStartupMessages(library(zoo))
invisible(Sys.setlocale("LC_TIME", "english"))
set.seed(1023)
mean.monthly.temp <- read.csv("mean-monthly-air-temperature-deg.csv", header=TRUE, stringsAsFactors = FALSE)
colnames(mean.monthly.temp) <- c("Date", "Temperature")
mmt <- mean.monthly.temp[complete.cases(mean.monthly.temp),]As you see above, I also cut out its last record as it reports just a description string without any temperature data associated.
Our dataset reports monthly data, so the frequency to be used in building the timeseries should be 12 (see ref. [3]). Anyway I want to verify that by means of the periodogram function made available by the TSA package. That could be useful in case you are uncertain about the frequency value to be specified.
temp.pdgram <- TSA::periodogram(mmt$Temperature, plot=TRUE)
(max.freq <- which.max(temp.pdgram$spec))## [1] 20(temp.frequency = 1/temp.pdgram$freq[max.freq])## [1] 12As we can see the frequency is confirmed to be equal to 12. I then create the ts time series object and plot it.
mean.monthly.ts <- ts(mmt$Temperature, start=c(1920,1), end=c(1939,12), frequency=temp.frequency)
tt <- time(mean.monthly.ts)
plot(mean.monthly.ts)
grid(col='blue',lty="dotted")
In my previous posts, I outlined the exploratory analysis of this dataset and also verified there are not any outliers or structural changes.
The following function collects summaries, basic metrics, stationarity and normality tests.
compute_tests <- function(timeserie) {
summary <- summary(timeserie)
sd <- sd(timeserie)
skew <- TSA::skewness(timeserie)
kurt <- TSA::kurtosis(timeserie)
adf.test <- tseries::adf.test(timeserie) # test for stationarity
jb.test <- tseries::jarque.bera.test(timeserie) # test for normality
lillie.test <- nortest::lillie.test(timeserie) # test for normality
list("summary" = summary,
"sd" = sd,
"skew" = skew,
"kurtosis" = kurt,
"adf" = adf.test,
"jarque" = jb.test,
"lillie" = lillie.test)
}Please see my previous post for comments about the compute_tests() results related the time series under analysis.
I build the Holt-Winters model with additive seasonal component. In doing that I take advantage of the forecast package hw() function.
mean.monthly.hw <- hw(mean.monthly.ts, seasonal = c("additive"))
fit.ts <- mean.monthly.hw$fitted
merge.ts.1 <- merge(as.zoo(fit.ts), as.zoo(mean.monthly.ts))
plot(merge.ts.1, screen=1, lty=c("dotted", "solid"), col= c("blue", "red"), ylim=c(30,75))
legend("topleft", c("fitted", "original"), lty=c("dotted", "solid"), col= c("blue", "red"), bty='n')
accuracy(mean.monthly.ts, fit.ts)## ME RMSE MAE MPE MAPE ACF1
## Test set 0.01354404 2.257418 1.773352 -0.002677307 3.787264 0.1643867
## Theil's U
## Test set 0.5486537The accuracy() function shows metrics to judge the evaluate the error between original time series and the fit just computed.
The hw() function returns also the forecast for next two periods that can be easily plot as herein shown.
plot(mean.monthly.hw)
The residual time series plots are herein shown.
resid.ts <- resid(mean.monthly.hw)
plot(resid.ts, type='l')
Residual analysis follows.
par(mfrow=c(1,3))
acf(resid.ts, lag.max=24)
pacf(resid.ts, lag.max=24)
McLeod.Li.test(y=resid.ts, plot=TRUE, main="McLeod-Li test")
LjungBoxTest(resid.ts, lag.max=24, k=14)## m Qm pvalue
## 1 1.79 1.809127e-01
## 2 2.14 1.432846e-01
## 3 5.38 2.034426e-02
## 4 6.04 1.400916e-02
## 5 7.16 7.446872e-03
## 6 7.23 7.164963e-03
## 7 7.23 7.163330e-03
## 8 7.25 7.100313e-03
## 9 11.45 7.168390e-04
## 10 13.76 2.077615e-04
## 11 23.24 1.430498e-06
## 12 23.81 1.064258e-06
## 13 24.95 5.893226e-07
## 14 24.96 5.854828e-07
## 15 25.69 4.007283e-07
## 16 26.36 1.890352e-06
## 17 26.93 6.078111e-06
## 18 27.75 1.398970e-05
## 19 27.76 4.056869e-05
## 20 27.81 1.019414e-04
## 21 29.73 1.063930e-04
## 22 29.75 2.335901e-04
## 23 35.93 4.085456e-05
## 24 37.91 3.938874e-05Some comments on aboveshown results.
based on the McLeod.Li test, the null hypothesis of homoscedasticity cannot be rejected.
the PACF plot indicates some autocorrelation is yet present at lags 10 and 12.
the LjungBox test confirms the residuals are not indipendentely distributed, i.e. they are correlated.
Furthermore, qq-plots and tests results are herein given.
qqnorm(resid.ts)
qqline(resid.ts)
resid.ts.test <- compute_tests(resid.ts)
resid.ts.test## $summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -7.03000 -1.33800 0.25280 -0.01354 1.42500 5.37400
##
## $sd
## [1] 2.262095
##
## $skew
## [1] -0.3083657
##
## $kurtosis
## [1] 0.1218341
##
## $adf
##
## Augmented Dickey-Fuller Test
##
## data: timeserie
## Dickey-Fuller = -5.9643, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $jarque
##
## Jarque Bera Test
##
## data: timeserie
## X-squared = 3.952, df = 2, p-value = 0.1386
##
##
## $lillie
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: timeserie
## D = 0.059534, p-value = 0.03856Based on Jarque-Bera test, we cannot reject the hypothesis of normality. On the contraty, based on Lillie test, we should reject at 5% p-value. However the Lillie test p-value is only slightly below 5%, a border-line scenario whereby it is hard to conclude.
I then take advantage of the forecast package auto.arima() function to investigate further the structure of the residual time series.
resid.auto.arima <- auto.arima(resid.ts, stepwise=FALSE, approximation=FALSE)
resid.auto.arima## Series: resid.ts
## ARIMA(1,0,0)(0,0,2)[12] with zero mean
##
## Coefficients:
## ar1 sma1 sma2
## 0.1836 -0.1547 0.1992
## s.e. 0.0640 0.0640 0.0731
##
## sigma^2 estimated as 4.704: log likelihood=-525.44
## AIC=1058.87 AICc=1059.04 BIC=1072.8An autoregressive model ARIMA(1,0,0)(0,0,2)[12] is the result of the auto.arima() procedure. The coefficients are significative.
par(mfrow=c(1,3))
acf(resid.auto.arima$resid, lag.max=24)
pacf(resid.auto.arima$resid, lag.max=24)
McLeod.Li.test(y=resid.auto.arima$resid, plot=TRUE, main="McLeod-Li test")
LjungBoxTest(resid.auto.arima$resid, lag.max=24, k=3)## m Qm pvalue
## 1 0.73 0.39413662
## 2 1.05 0.30607153
## 3 3.62 0.05692578
## 4 3.63 0.05681420
## 5 4.88 0.08731751
## 6 5.05 0.16824165
## 7 5.07 0.28005853
## 8 5.32 0.37808812
## 9 7.65 0.26519691
## 10 8.36 0.30223484
## 11 8.41 0.39485255
## 12 9.63 0.38097536
## 13 11.98 0.28607176
## 14 12.71 0.31251168
## 15 14.58 0.26524155
## 16 14.99 0.30773590
## 17 15.83 0.32415621
## 18 17.65 0.28177094
## 19 17.75 0.33852073
## 20 17.81 0.40066841
## 21 18.05 0.45237073
## 22 19.10 0.45065724
## 23 19.10 0.51503125
## 24 19.91 0.52708603compute_tests(resid.auto.arima$resid)## $summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -6.050000 -1.360000 0.117400 -0.006751 1.258000 5.332000
##
## $sd
## [1] 2.159659
##
## $skew
## [1] -0.09799722
##
## $kurtosis
## [1] -0.1884668
##
## $adf
##
## Augmented Dickey-Fuller Test
##
## data: timeserie
## Dickey-Fuller = -6.2012, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $jarque
##
## Jarque Bera Test
##
## data: timeserie
## X-squared = 0.73934, df = 2, p-value = 0.691
##
##
## $lillie
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: timeserie
## D = 0.039789, p-value = 0.4682Some comments on aboveshown results.
the residuals are uncorrelated
we cannot reject the null hypothesis on residuals normal distribution
we cannot reject the null hypothesis of homoscedasticity on residuals
Our time serie can be so represented:
\[ \begin{equation} \begin{cases} y_{t}\ =\ HoltWinters(x_{t}) +\ z_{t} \\ \\ z_{t}\ =\ 0.184 z_{t-1}\ +\ \eta_{t} \\ \\ \eta_{t} =\ -0.155 \eta_{t-12}\ +\ 0.199 \eta_{t-24}\ +\ \epsilon_{t} \\ \\ \epsilon_{t} \sim WN(0,\ 4.704) \\ \end{cases} \end{equation} \]
I then simulate such model and adding the resulting time series to the fitted one.
resid.auto.sim <- simulate(resid.auto.arima, future=FALSE)
resid.auto.sim.ts <- ts(coredata(resid.auto.sim), start=c(1920,1), end=c(1939,12), frequency=temp.frequency)
fit.ts.2 <- fit.ts + resid.auto.sim.ts
accuracy(mean.monthly.ts, fit.ts.2)## ME RMSE MAE MPE MAPE ACF1
## Test set 0.04880902 3.305397 2.665326 -0.1431367 5.653688 0.1837114
## Theil's U
## Test set 0.6653601merge.ts.2 <- merge(as.zoo(fit.ts.2), as.zoo(mean.monthly.ts))
plot(merge.ts.2, screen=1, lty=c("dotted", "solid"),col= c("blue", "red"), ylim=c(30,75))
legend("topleft", c("fitted", "original"), lty=c("dotted", "solid"), col= c("blue", "red"), bty='n')
Above accuracy() function metrics reveal that the second model remainder time series (which includes the ARIMA model based simulated component) shows a lower goodness of fit with respect the first fit I computed. However it is the result of a modeling which better capture the structure of the original time series as its residuals have been verified to be white noise.
Ultimately, some analysis on the remainder associated to the second fit which is based on the sum up of the Holt-Winters fit and its residuals ARIMA model simulation
resid.fit.2 <- mean.monthly.ts - fit.ts.2
par(mfrow=c(1,3))
acf(resid.fit.2, lag.max=24)
pacf(resid.fit.2, lag.max=24)
(m <- McLeod.Li.test(y=resid.fit.2, plot=TRUE, main="McLeod-Li test"))
## $p.values
## [1] 0.04006783 0.08898026 0.18395676 0.29436724 0.36035431 0.47467339
## [7] 0.51894875 0.58279336 0.67602085 0.75908902 0.27755774 0.35141751
## [13] 0.42261536 0.49153540 0.56502060 0.62984926 0.69644707 0.75529170
## [19] 0.80694592 0.84348882 0.86073424 0.80100460 0.83261598LjungBoxTest(resid.fit.2, lag.max=24, k=0)## m Qm pvalue
## 1 0.16 0.6931029791
## 2 4.26 0.1188429403
## 3 4.57 0.2060677193
## 4 4.58 0.3335353928
## 5 4.58 0.4688092357
## 6 4.62 0.5933000926
## 7 7.78 0.3521959286
## 8 7.96 0.4378396140
## 9 8.50 0.4850467219
## 10 18.11 0.0531315869
## 11 27.46 0.0039096559
## 12 28.14 0.0052693065
## 13 32.55 0.0019885530
## 14 32.69 0.0031883769
## 15 32.82 0.0049761014
## 16 33.31 0.0067236693
## 17 34.81 0.0065936778
## 18 34.93 0.0096586581
## 19 34.93 0.0142357372
## 20 35.01 0.0200717009
## 21 36.45 0.0194762557
## 22 38.00 0.0183361097
## 23 50.43 0.0008096733
## 24 50.60 0.0011873245compute_tests(resid.fit.2)## $summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -8.52400 -2.22400 0.39180 -0.04881 2.27500 8.22800
##
## $sd
## [1] 3.311944
##
## $skew
## [1] -0.1947143
##
## $kurtosis
## [1] -0.338377
##
## $adf
##
## Augmented Dickey-Fuller Test
##
## data: timeserie
## Dickey-Fuller = -6.0738, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
##
##
## $jarque
##
## Jarque Bera Test
##
## data: timeserie
## X-squared = 2.6615, df = 2, p-value = 0.2643
##
##
## $lillie
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: timeserie
## D = 0.058825, p-value = 0.04302Some comments on aboveshown results.
the residuals shows correlations at lags {11, 12, 24}
based on Jarque-Bera test, we cannot reject the null hypothesis on residuals normal distribution
based on Lillie test, we can reject the null hypothesis on residuals normal distribution
based on McLeod-Li test, we reject the null hypothesis of homoscedasticity on residuals at lag 1. We cannot do the same for higher lag than one.
Anyway, since the second model residuals \(\epsilon_{t}\) are white noise, those remainder auto-correlation and heteroscedasticity appear to be more likely a specific consequence of the ARIMA simulation rather than a modeling issue.
Conclusions
When the time series shows evident seasonality, to build a Holt-Winter model is suggestable for a good fitting of the time series under analysis. I also introduced how easily forecasts could be obtained from.
Analysis of residuals has spot further structure in the Holt-Winters residuals whose ARIMA structure has been determined. The final residuals were verified to be white noise.
Two models remainders have been compaired against several metrics.
References
Note
Original post published on March 1st 2015 has been updated by substantial reworking and enhancing the analysis.
