Abstract

In this report we investigate a time-series dataset involving the monthly average minimum temperatures for the city of Melbourne, Australia, from the years 1981 to 1990. In particular, we look to address how to accurately model the dataset and to use this model to make predictions for the following year. Our investigation includes model building techniques with data transformation using the Box-Cox method, model selection among candidates based on the Akaike Information criterion and various diagnostic tests, and finally forecasting of future data points using our model. Based on our final tests we conclude that the final model we develop successfully describes the dataset and can be used to accurately predict future data points.

Introduction

The dataset we are working with consists of the monthly average minimum temperatures for the city of Melbourne, Australia, provided by the Australian Bureau of Meteorology. It is a dataset spanning 10 years from 1981 to 1990, totaling 120 entries. Exploration of this dataset will provide an intuitive way of analyzing temperature patterns for the city. Our goal of this report is to be able to forecast the next 12 values of the dataset; that is, monthly averages for the following year.

Our model is constructed using RStudio software. In building our model we use various built-in techniques such as data transformation to make the dataset stationary, model selection through ACF/PACF analysis, and model diagnostic tests by performing residual analyses. The resulting model is one that has passed all of these tests and shows accurate forecasting ability in correspondence with our test dataset.

The Dataset

Let us first take a look at the plot of the dataset:

temp.csv<-read.table("monthly-avg-min-temperatures.csv", sep=',', header=TRUE)
temp_ts<-ts(temp.csv['X0'], start=1981, frequency=12)
ts.plot(temp_ts, main="Monthly Average Minimum Temperatures",
        ylab="Temperature (C)", xlab="Year")

We observe that the dataset is clearly non-stationary due to its pronounced seasonality. We also notice that there seems to be slight changes in variance throughout the data. There does not seem to be any apparent trend or other odd behaviors.

Since we want to build up our model, we will be truncating this data by the last 12 entries, leaving those entries to be used for model diagnostics later on. Let us call this new dataset \(Y_t,\;t=1,2,...,108\).

temp.train<-temp_ts[c(1:108)]  # training dataset
temp.test<-temp_ts[c(109:120)]  # test dataset
ts.plot(temp.train, main="Plot of Truncated Data Y_t")

acf(temp.train, main="ACF of Y_t")

The plot of the dataset and of its ACFs are a big indicator for adjusting the dataset to make it stationary. Since variance is uneven, we may try to perform a variance stabilizing transformation on the data. The ACFs are large and periodic, so we may try differencing on the data.

Data Transformation

We first attempt to stabilize the variance of the data by performing a Box-Cox transformation:

library(MASS)
bcTransform<-boxcox(temp.train~as.numeric(1:length(temp.train)))

lambda<-bcTransform$x[which(bcTransform$y==max(bcTransform$y))]
lambda
## [1] 0.7070707

Using the bcTransform() command gives the value \(\lambda=0.7070707\). Now we transform the data using the equation \(U_t=\frac{1}{\lambda}(Y_t^{\lambda}-1)\), where \(U_t\) is our new transformed data.

To see if variance has been stabilized, let us compare the plots of \(Y_t\) and \(U_t\):

temp.bc<-(1/lambda)*((temp.train^lambda)-1)
plot.ts(temp.train, main="Y_t")

plot.ts(temp.bc, main="U_t")

It is slight, but variance does look to be more stable after transformation, so we will stick with \(U_t\).

Next, we move on to differencing the data to remove the seasonality component. It is a monthly time-series dataset, so we will first difference the data at lag 12 and check variance to see if it is lower. Since trend is not apparent, we may not need to difference any more, in case we over-difference the data. Below is a plot of each differencing step along with the corresponding variance value below it:

plot.ts(temp.bc, main="U_t")

var(temp.bc)
## [1] 2.475097
temp.d<-diff(temp.bc, lag=12)
plot.ts(temp.d, main="U_t differenced at lag 12")

var(temp.d)
## [1] 0.4511216
temp.dd<-diff(temp.d, lag=1)
plot.ts(temp.dd, main="U_t differenced at lag 12 and lag 1")

var(temp.dd)
## [1] 0.7514186

We can see that after differencing once at lag 12, seasonality is no longer present and variance is lower. However, differencing once more at lag 1 increases variance, meaning we have over-differenced. Thus, we stick to differencing once at lag 12. Next we compare the ACFs of \(U_t\) and \(U_t\) differenced at lag 12:

acf(temp.bc, lag.max=50, main="ACF of U_t")

acf(temp.d, lag.max=50, main="ACF of U_t differenced at lag 12")

After differencing at lag 12, the ACFs show decay, indicating a stationary process. We can now use this to move on to model building for our data.

Model Building

In order to identify the possible parameters for our model, we will need ot analyze the ACF and PACF of the transformed and differenced data:

acf(temp.d, lag.max=50, main="ACF of U_t differenced at lag 12")

ACFs outside of confidence interval: lags 2, 12

pacf(temp.d, lag.max=50, main="PACF of U_t differenced at lag 12")

PACFs outside of confidence interval: lags 2, 12

Then we have our list of candidate models to try:

Model Selection

We now use the arima() function to obtain coefficient values for the models, as well as use the AICc() function to compare AICc between the models:

ar1<-arima(temp.bc, order=c(0,0,2), seasonal=list(order=c(0,1,1), period=12),
           method="ML")
ar2<-arima(temp.bc, order=c(2,0,0), seasonal=list(order=c(1,1,0), period=12),
           method="ML")
ar3<-arima(temp.bc, order=c(2,0,2), seasonal=list(order=c(1,1,1), period=12),
           method="ML")
library("MuMIn")
## Warning: package 'MuMIn' was built under R version 4.0.5
ar1$coef  # SMA: s=12, d=0, D=1, q=2, Q=1
##        ma1        ma2       sma1 
##  0.1832781  0.2628922 -0.8757161
AICc(ar1)
## [1] 160.4064
ar2$coef  # SAR: s=12, d=0, D=1, p=2, P=1
##        ar1        ar2       sar1 
##  0.1280750  0.2767044 -0.4784034
AICc(ar2)
## [1] 170.8895
ar3$coef  # SARIMA: s=12, d=0, D=1, p=q=2, P=Q=1
##         ar1         ar2         ma1         ma2        sar1        sma1 
##  0.48605643 -0.10250108 -0.28828833  0.31863585  0.08066783 -0.99987692
AICc(ar3)
## [1] 165.4036

We want our model to have a low AICc; we can see that the SMA and SARIMA models have the two lowest AICcs, so we choose these two as our candidate models to perform further model diagnostics on.

Then our candidate models are

Model Diagnostics for Model A

We now take a look at the residuals of model A, checking for normality assumption:

fitA<-arima(temp.bc, order=c(2,0,2), seasonal=list(order=c(1,1,1), period=12),
           method="ML")
resA<-residuals(fitA)
hist(resA, col="red", density=30, prob=TRUE, main="Histogram of Model A Residuals")
m<-mean(resA)
std<-sqrt(var(resA))
curve(dnorm(x,m,std), add=TRUE)

plot.ts(resA, main="Plot of Model A Residuals")
fitA2<-lm(resA~as.numeric(1:length(resA))); abline(fitA2, col="red")
abline(h=mean(resA), col="blue")

qqnorm(resA, main="Normal Q-Q Plot for Model A")
qqline(resA, col="blue")

There does not seem to be anything out of the ordinary with the residual plots; no seasonality, no trend, and variance is stable. Both histogram and normal Q-Q plots look okay as well.

acf(resA, lag.max=50, main="ACF of Model A Residuals")

pacf(resA, lag.max=50, main="ACF of Model A Residuals")

ACFs and PACFs are all within confidence intervals.

shapiro.test(resA)
## 
##  Shapiro-Wilk normality test
## 
## data:  resA
## W = 0.98283, p-value = 0.1786
Box.test(resA, lag=11, type=c("Box-Pierce"), fitdf=2)
## 
##  Box-Pierce test
## 
## data:  resA
## X-squared = 5.7693, df = 9, p-value = 0.7628
Box.test(resA, lag=11, type=c("Ljung-Box"), fitdf=2)
## 
##  Box-Ljung test
## 
## data:  resA
## X-squared = 6.3756, df = 9, p-value = 0.7018
Box.test(resA^2, lag=11, type=c("Ljung-Box"), fitdf=0)
## 
##  Box-Ljung test
## 
## data:  resA^2
## X-squared = 14.514, df = 11, p-value = 0.2058

All p-values are greater than 0.05; model A has passed all diagnostic checks.

Model Diagnostics for Model B

We now do the same checks for Model B:

fitB<-arima(temp.bc, order=c(0,0,2), seasonal=list(order=c(0,1,1), period=12),
           method="ML")
resB<-residuals(fitB)
hist(resB, col="green", density=30, prob=TRUE, main="Histogram of Model B Residuals")
m<-mean(resB)
std<-sqrt(var(resB))
curve(dnorm(x,m,std), add=TRUE)

plot.ts(resB, main="Plot of Model B Residuals")
fitB2<-lm(resB~as.numeric(1:length(resB))); abline(fitB2, col="red")
abline(h=mean(resB), col="blue")

qqnorm(resB, main="Normal Q-Q Plot for Model B")
qqline(resB, col="blue")

Residual plots look okay as there is no seasonality, no trend, and variance is stable. Both histogram and normal Q-Q plots look fine.

acf(resB, lag.max=50, main="ACF of Model B Residuals")

pacf(resB, lag.max=50, main="ACF of Model B Residuals")

ACFs and PACFs are within confidence intervals.

shapiro.test(resB)
## 
##  Shapiro-Wilk normality test
## 
## data:  resB
## W = 0.98086, p-value = 0.1225
Box.test(resB, lag=11, type=c("Box-Pierce"), fitdf=3)
## 
##  Box-Pierce test
## 
## data:  resB
## X-squared = 5.7438, df = 8, p-value = 0.6759
Box.test(resB, lag=11, type=c("Ljung-Box"), fitdf=3)
## 
##  Box-Ljung test
## 
## data:  resB
## X-squared = 6.2709, df = 8, p-value = 0.6169
Box.test(resB^2, lag=11, type=c("Ljung-Box"), fitdf=0)
## 
##  Box-Ljung test
## 
## data:  resB^2
## X-squared = 14.628, df = 11, p-value = 0.2002

All p-values are greater than 0.05; Model B has also passed all diagnostic checks.

Both of our models have shown to pass the diagnostic tests. To choose between the two, we look at which model is parsimonious; that is, the model that requires the least amount of parameters to estimate. Thus we will be choosing Model B over model A for forecasting of the data.

Our final model for \(U_t\) follows \(\text{SARIMA}(0,0,2)\times(0,1,1)_{12}\\\nabla_{12}U_t=(1+0.183B+0.263B^2)(1-0.876B^{12})Z_t\)

Forecasting

We now turn to forecasting with our model B using predict():

pred<-predict(fitB, n.ahead=12)
U<-pred$pred + 2*pred$se
L<-pred$pred - 2*pred$se
ts.plot(temp.bc, xlim=c(1, length(temp.bc)+12), ylim=c(min(temp.bc), max(U)),
        main="Forecast of U_t using Model B")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
points((length(temp.bc)+1):(length(temp.bc)+12), pred$pred, col="red")

pred.orig<-exp(log(lambda*pred$pred+1)/lambda)
U<-exp(log(lambda*U+1)/lambda)
L<-exp(log(lambda*L+1)/lambda)
ts.plot(temp.train, xlim=c(1, length(temp.train)+12), ylim=c(min(temp.train), max(U)),
        main="Forecast of Y_t using Model B")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
points((length(temp.train)+1):(length(temp.train)+12), pred.orig, col="red")

ts.plot(temp.train, xlim=c(60, length(temp.train)+12), ylim=c(5, max(U)),
        main="Zoomed in forecast of Y_t")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
points((length(temp.train)+1):(length(temp.train)+12), pred.orig, col="red")

temp_ts<-ts(temp.csv['X0'])
ts.plot(temp_ts, xlim=c(60, length(temp.train)+12), ylim=c(5, max(U)),
        main="Forecast of Y_t including test set")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
points((length(temp.train)+1):(length(temp.train)+12), pred.orig, col="red")

We can see that when plotting the predicted values on top of the test set, the test set is within the prediction intervals. Our model looks to be fairly accurate.

Conclusion

To conclude, we were able to achieve our goal of forecasting future temperatures using this dataset. Our model was able to pass all diagnostic checks and produce an accurate forecast. When testing the model on our test set, we were able to see that the test set data entries fell within the produced prediction intervals. Again the model for our data is \(\nabla_{12}U_t=(1+0.183B+0.263B^2)(1-0.876B^{12})Z_t\).

References

“Australian Government Bureau of Meteorology.” Australia’s Official Weather Forecasts & Weather Radar - Bureau of Meteorology.