TIMES SERIES PROJECT: MELBOURNE TEMPERATURES
ThÃffÃ???’Ãf,Ã,©o Dufort S3593458, Simon Terry S3610872, Siddharth Jain S3616023
26/05/2017
INTRODUCTION
In this report we will perform a time series analysis of Melbourne’s average monthly temperatures. Firstly we will examine the time series and attempt to specify a model to fit historic temperature variations. After this we will use the selected model to forecaste the next 10 years of Melbourne’s monthly average temperatures.
Note: For the purposes of assessment, all chunks of R code will be visible.
DATA DESCRIPTION
Data set is Melbourne’s monthly average temperatures from January 1841 to September 2013.
Source : http://berkeleyearth.lbl.gov/auto/Local/TAVG/Text/37.78S-144.41E-TAVG-Trend.txt
% This file contains an extracted local summary of land-surface % temperature results produced by the Berkeley Earth averaging % method for the location: % % 37.78 S, 144.41 E % % The Berkeley Earth method takes temperature observations from a large % collection of weather monitoring stations and produces an estimate of % the underlying global temperature field across all of the Earth’s % land areas. Once this temperature field has been generated, it is % possible to estimate the temperature evolution of individual locations % simply by sampling the field at the locaiton in question. This % file contains such a local estimate. % % % The current dataset presented here is described as: % % Berkeley Earth analysis for mean temperature on complete dataset % % % This analysis was run on 12-Oct-2013 00:45:15 % % Global results are based on 39348 time series % with 15244148 data points % % The current location is characterized by: % Country: Australia % Nearby Cities: Melbourne, Geelong, Ballarat, Cranbourne, Melton, Sunbury % Percent water in local neighborhood: 16.2 % % Temperature stations within 200 km: 145 % Temperature obeservations within 200 km: 45302 %
METHODOLOGY
In order to specify the model we first need to examine the data for stationarity. Obviously, monthly average temperatures are going to exhibit seasonality. The first step is to seasonally adjust the time series.
After seasonality has been acccounted for we need to check the resulting seasonally adjusted series for autocorrelation, trend and stationarity. If necessary we will difference the data series.
The next step will be to create candidate models to account for the remaining autocorrelation. Candidate models with remaining autocorrelation will be elimianated. Remaining models will be subjected to coefficient testing and models with insignificant coefficients will be removed. The candidate models that are still considered viable will be compared with residual analysis and the best model selected as our final model.
The final model will be used to predict the next 120 months (10 years) of Melbourne’s average monthly temperatures.
TIME SERIES ANALYSIS
Importing
Importing the libraries
The packages needed for analysis are loaded in R.
rm(list=ls())
library(TSA) # To have nice ACFs
library(lmtest) #For coeftest()
library(fUnitRoots) #For unit root test
library(readr) #To import the data set
library(FitAR) #For residual AnalysisImporting the data
The Melbourne temperatures’ data set is imported.
BerkleyTempNewDataset <- read_csv("C:/Users/siddh/Desktop/Study/Sem_2/Time Series Analysis/Project/BerkleyTempNewDataset.csv",
col_types = cols(`Anomaly,` = col_skip(),`Anomaly,_1` = col_skip(),
`Anomaly,_2` = col_skip(), `Anomaly,_3` = col_skip(),
`Anomaly,_4` = col_skip(), Unc. = col_skip(),
`Unc.,` = col_skip(), `Unc.,_1` = col_skip(),
`Unc.,_2` = col_skip(), `Unc.,_3` = col_skip(),
X3 = col_skip()))
Temperature<-ts(BerkleyTempNewDataset$X5,start = 1841,frequency = 12)
class(Temperature)## [1] "ts"Temperature is now a times series, with one temperature value showing average temperature per month from January of 1841 to September of 2013.
In this analysis, the focus is made on the 1950-2013 period. We selected this time pan as it exhibited an upward trend not apparent in earlier years. Consequently earlier observations are taken out of the data set. Here is a plot of the new time series:
Temperatures <- window(Temperature , 1950)
plot(Temperatures, ylab='Temperature (°C)',main='Melbourne Temperatures between 1950 and 2013 ')
Month=c('J','F','M','A','M','J','J','A','S','O','N','D')
points(Temperatures,pch=Month) 
At this point, seasonality is very obvious, with peaks in temperature for the months of January and February, and lowest values for the months of July and August. This comes as no surprise as Melbourne has 4 seasons (in a day according to legend), with warmer weather in summer and colder weather in winter. However this seasonality makes it hard to spot any other kind of trend.
Seasonal Component
A look at the Auto-Correlation Function (ACF) and Partial Auto-Correlation Function (PACF):
acf(as.vector(Temperatures),lag.max=36, ci.type="ma", main='Temperatures') 
acf(as.vector(Temperatures),lag.max=36,type="partial",ci.type="ma",main='Temperatures',ylab="PACF")
The ACF confirms the obvious, with very high positive correlations with lags 11, 12, 13 and high negative correlations with lags 5,6,7 etc… This clearly indicates a yearly seasonal pattern. This pattern fails to decrease over time, which indicates that a seasonal differencing using lag 12 is in order. At this point, it is becoming clear that a SARIMA model has to be used for modelling.
The PACF shows high correlations on the first lags, with a geometric decay that once again alternates between positive and negative.
As explained, a differencing is done with lag 12, and the results are then plotted:
plot(diff(Temperatures,lag=12),ylab='First Difference of Melbourne Temperature (°C)',xlab='Years',
main = "First Differences of Melbourne Temperatures Over Time") #The first 12 values are lost.
At first glance, the seasonality seems to have disappeared. however, when adding the months’ names on each observation…
plot(diff(Temperatures,lag=12),ylab='First Difference of Melbourne Temperature (°C)',xlab='Years',
main = "First Differences of Melbourne Temperatures Over Time") #The first 12 values are lost.
points(diff(Temperatures),pch=Month) 
It now seems that October, November and December often have higher values than average, and the opposite is observed for February, March, April and May. It indicates that a certain degree of seasonality is still present in the data.
the corresponding ACF and PACF confirm this finding.
par(mfrow=c(1,2))
acf(as.vector(diff(Temperatures,lag=12)),lag.max=36,ci.type="ma",
main="Sample ACF of First Differences
of Melbourne Temperature")
acf(as.vector(diff(Temperatures,lag=12)),lag.max=36,type="partial",ci.type="ma",
main="Sample PACF of First Differences
of Melbourne Temperature",ylab="PACF")
Indeed, the ACF show a very strong correlation with lag 12, and the PACF with lag 12, 24, 36 etc.. this is typical of a Seasonal Moving Average (SMA).
Hence, a SARIMA(0,0,0)x(0,1,1)_12 model is built from the raw data, and a plot of the residuals is generated.
#SARIMA(0,0,0)X(0,1,1)_12
par(mfrow=c(1,1))
mx.Temperatures = arima(Temperatures,order=c(0,0,0),seasonal=list(order=c(0,1,1), period=12))
res.mx = residuals(mx.Temperatures)
plot(res.mx,xlab='Time',ylab='Residuals',main="Time series")
abline(h=0,col="blue",lwd=2)
The plot shows no seasonal pattern, however the residuals seem to follow each other, which is a sign of autocorrelation. Finally, they gravitate around the zero line, which is a stationary behaviour.
The corresponding ACF and PACF are plotted:
par(mfrow=c(1,2))
acf(as.vector(res.mx),lag.max=36,ci.type="ma", main="Sample ACF of the Residuals")
acf(as.vector(res.mx),lag.max=36,ci.type="ma",type="partial",main="Sample PACF of the Residuals",ylab="PACF")
Finally, it appears that no seasonality is left . Indeed, neither the ACF nor PACF show significant correlation with lag 12. The focus can now be made on the ARIMA component of the SARIMA model.
ARIMA Component: Finding Candidates
The ACF showed significant correlations with lags 1 & 2 (potentially lag 3 as well, although it might simply be a geometric decay decreasing very quickly) and the PACF showed significant correlation with lag 1, maybe lag 2. There is no significance in late lags and barely any geometric pattern, hence no sign of non-stationarity whatsoever. To confirm, a Unit Root Test is run:
order = ar(diff(res.mx))$order
adfTest(res.mx, lags = order, title = NULL,description = NULL)@test$p.value##
## 0.01P value is 0.01 so the non-stationarity hypothesis can be rejected at the 5% level. Consequently, it is decided to model the following list of candidates:
*M1 = SARIMA(0,0,1)x(0,1,1)_12 *M2 = SARIMA(1,0,0)x(0,1,1)_12 *M3 = SARIMA(1,0,1)x(0,1,1)_12 *M4 = SARIMA(0,0,2)x(0,1,1)_12 *M5 = SARIMA(2,0,0)x(0,1,1)_12 *M6 = SARIMA(2,0,1)x(0,1,1)_12 *M7 = SARIMA(1,0,2)x(0,1,1)_12 *M8 = SARIMA(2,0,2)x(0,1,1)_12
M1 : SARIMA(0,0,1)X(0,1,1)_12
m1.Temperatures = arima(Temperatures,order=c(0,0,1),seasonal=list(order=c(0,1,1),
period=12),method="ML")
m1.Temperatures_css = arima(Temperatures,order=c(0,0,1),seasonal=list(order=c(0,1,1),
period=12),method="CSS")
res.m1 = residuals(m1.Temperatures)
par(mfrow=c(1,2))
acf(as.vector(res.m1),lag.max=36,ci.type="ma", main="Sample ACF of the Residuals")
acf(as.vector(res.m1),lag.max=36,ci.type="ma",type="partial",main="Sample PACF of the Residuals",ylab="PACF")
Both ACF and PACF show a significant correlation with lag 2. Thus, this model is dropped from the list of candidates.
M2 : SARIMA(1,0,0)X(0,1,1)_12 
Some lags (eg. lag 9 in ACF, lag 3 and 9 in PACF) are still significant, although barely. This is a satisfactory model, however it is easily improved by adding one more term to the ARMA component (as seen next).
M3 : SARIMA(1,0,1)X(0,1,1)_12 
Once more, lag 9 is barely significant on both the ACF and PACF (and the presumed correlation is under 0.10), however this time the other lags are very clearly insignificant. This model is the best one so far.
See below the coefficients computed for this model, using the Maximum Likelihood method (ML) and Cumulative Sum of Square method (CSS) respectively.
coeftest(m3.Temperatures)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.698947 0.092926 7.5215 5.414e-14 ***
## ma1 -0.489797 0.111701 -4.3849 1.160e-05 ***
## sma1 -0.935175 0.018245 -51.2567 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(m3.Temperatures_css)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.546591 0.076759 7.1209 1.073e-12 ***
## ma1 -0.303075 0.082083 -3.6923 0.0002222 ***
## sma1 -0.875868 0.019140 -45.7615 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1For both methods, all coefficients are significant at the 1% level (This validates the decision to drop models M1 and M2). For both results, the AR component is bigger than 0.5, the MA one is between -0.30 and -0.5 and the SMA one is close to -0.90.
M4 : SARIMA(0,0,2)X(0,1,1)_12 
This time, the plots both show significant correlation with lag 3. It appears models without AR() component do not fit the data properly. This model is excluded from the set of candidates.
M5 : SARIMA(2,0,0)X(0,1,1)_12 
Results are similar to the findings for M3. This is a good model.
See below the coefficients computed for this model, using the Maximum Likelihood ML and CSS respectively.
coeftest(m5.Temperatures)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.219672 0.036641 5.9952 2.033e-09 ***
## ar2 0.108330 0.036675 2.9538 0.003139 **
## sma1 -0.927578 0.017545 -52.8700 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(m5.Temperatures_css)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.238073 0.036820 6.4659 1.007e-10 ***
## ar2 0.115661 0.036530 3.1662 0.001545 **
## sma1 -0.876385 0.019121 -45.8327 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Once again, all coefficients are significant at the 1% level for both methods, however this time the value computed for each coefficient is very similar from one method to the other. Hence, the first AR component, is slightly above 0.2, the second is slightly above 0.10, and the SMA coefficient is close to -0.90.
M6 : SARIMA(2,0,1)X(0,1,1)_12 
The ACF and PACF for this model are extremely similar to the previous ones, except that correlation to lag 9 is slightly smaller.
Looking at the coefficient computed via ML and CSS.
coeftest(m6.Temperatures)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.790794 0.287010 2.7553 0.005864 **
## ar2 -0.032785 0.100232 -0.3271 0.743602
## ma1 -0.577788 0.281807 -2.0503 0.040336 *
## sma1 -0.936752 0.019255 -48.6496 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(m6.Temperatures_css)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.589414 0.186138 3.1665 0.001543 **
## ar2 0.025110 0.067385 0.3726 0.709421
## ma1 -0.357941 0.181466 -1.9725 0.048552 *
## sma1 -0.878899 0.018976 -46.3162 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1This time, the second AR coefficient is not significant with a P-value > 0.70 in both cases. Consequently the model is excluded from the list of candidate models.
M7 : SARIMA(1,0,2)X(0,1,1)_12 
This model gives the same ACF and PACF as model M6 seen previously.
Looking at the computed coefficients:
coeftest(m7.Temperatures)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.740878 0.149921 4.9418 7.741e-07 ***
## ma1 -0.528040 0.153733 -3.4348 0.000593 ***
## ma2 -0.021525 0.065014 -0.3311 0.740581
## sma1 -0.936553 0.018947 -49.4308 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(m7.Temperatures_css)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.511635 0.087987 5.8149 6.068e-09 ***
## ma1 -0.284463 0.091004 -3.1258 0.001773 **
## ma2 0.045501 0.040080 1.1353 0.256269
## sma1 -0.876334 0.019088 -45.9097 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Once again, not all components are significant: The second MA coefficient is insignificant in both cases with P-value > 0.25. This model is excluded form the list of candidates. As ARMA(2,1) and ARMA(1,2) have insignificant coeffecients, ARMA (2,2) is an overfitted on model, hence they are excluded from the candidate models.
To summarise, the best two models at this stage of the analysis are: *M3: SARIMA(1,0,1)X(0,1,1)_12 *M5: SARIMA(2,0,0)X(0,1,1)_12
They will both be run through residual analysis.
Residual Analysis
Find below the function used for residual analysis. It displays a plot of the standardised residuals, with the computed proportion of outliers (more than 2 units away from zero.). It also shows an histogram of their distribution, along with a QQ plot and the computed p-value for the Shapiro test. Finally the p-value of the Ljung-Box Test is plotted for the first 40 lags. The ACF and PACF of the standardised residuals are omitted as they are identical to those of the regular residuals analysed in the previous section.
residual.analysis <- function(model, std = TRUE){
if (std == TRUE){
res.model = rstandard(model)
}else{
res.model = residuals(model)
}
par(mfrow=c(2,2))
plot(res.model,type='o',ylab='Standardised Residuals', main="Time Series Plot of Std Residuals")
abline(h=0,col="blue",lwd=2)
abline(2,0,col="green",lty=2)
abline(-2,0,col="green",lty=2)
abline(-3,0,col="red",lwd=2)
abline(3,0,col="red",lwd=2)
prop.outliers=((sum(res.model>2)+sum(res.model<(-2)))/ sum(res.model<30))
print(merge("Proportion of Residuals Greater Than 2 Standard Deviations",prop.outliers))
hist(res.model,main="Histogram of Standardised Residuals")
qqnorm(res.model,main="QQ plot of Standardised Residuals")
qqline(res.model, col = 2)
print(shapiro.test(res.model))
k=0
LBQPlot(res.model, lag.max = 40, StartLag = k + 1, k = 0, SquaredQ = FALSE)
}M3 : SARIMA(1,0,1)X(0,1,1)_12
Pushing the M3 model into the residual.analysis() function gives the following results.
residual.analysis(m3.Temperatures)## x y
## 1 Proportion of Residuals Greater Than 2 Standard Deviations 0.05490196##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.993, p-value = 0.001188## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple
## of shorter object length
The standardised residuals seem to behave like white noise. A few observations are greater than 3, however the overall proportion of residuals more than 2 standard deviations away from zero is 0.055, extremely close to the expected 0.05. Moreover, the histogram shows an almost symetrical distribution around 0, and tails off on both ends, which once again is typical behaviour of white noise. The QQ plot only confirms this, with sychronised sample and theorical quantiles, except at the upper end where the residuals appear slightly above the QQ line. Despite these good looking graphs, the Shapiro test of normality can be rejected at the 5% level with a p-value of 0.0012. Finally, the analysis ends on a high note with all p-values from the Ljung-Box Test clearly above 0.4, far away from the 0.05 line. The residuals from the SARIMA(1,0,1)x(0,1,1)_12 are overall very satisfactory.
M5 : SARIMA(2,0,0)X(0,1,1)_12
Below are the results from the standard residuals analysis on the M5 model.
residual.analysis(m5.Temperatures)## x y
## 1 Proportion of Residuals Greater Than 2 Standard Deviations 0.05751634##
## Shapiro-Wilk normality test
##
## data: res.model
## W = 0.99269, p-value = 0.0008276## Warning in (ra^2)/(n - (1:lag.max)): longer object length is not a multiple
## of shorter object length
The results are extremely similar to the findings from the M3 (SARIMA(1,0,1)X(0,1,1)_12) analysis. The plot seems almost identical, athlough the proportion of residuals more than 2 standard deviations away from zero is slightly larger (0.058). The histogram and QQ plot are also very reminiscent of M3. The Shapiro test rejects normality as well, although with an even smaller p-value (0.0008). Finally, the p-values from the Ljung-Box Test draw almost the same pattern, only with some smaller p-values. Indeed, the tests for lag 1 and 2 result in p-values clearly under 0.4. To conclude, this model gives excellent results as well, however not as good as M3.
AIC/BIC
Finally, another method to determine which model is the most appropriate is to look at the AIC and BIC scores. Reminder: The smaller the score, the better the model.
AIC(m3.Temperatures,m5.Temperatures)## df AIC
## m3.Temperatures 4 2157.015
## m5.Temperatures 4 2160.631The AIC gives a slightly smaller score for the M3 (SARIMA(1,0,1)x(0,1,1)_12) model with 2157 against 2161 for M5 (SARIMA(2,0,0)X(0,1,1)_12).
BIC(m3.Temperatures,m5.Temperatures)## df BIC
## m3.Temperatures 4 2175.511
## m5.Temperatures 4 2179.127The BIC seems to confirm it with 2176 against 2179 for M3 (SARIMA(1,0,1)x(0,1,1)_12) and and M5 (SARIMA(2,0,0)X(0,1,1)_12) respectively.
Overall, even though both models seems to perform very well, the study shows that the SARIMA(1,0,1)x(0,1,1)_12 is the most fitted for Melbourne’s temperatures data. Therfore it is selected for predictions.
PREDICTIONS
Find below a plot of the monthly temperatures in Melbourne from January 2000 to September 2013, followed by ten years of predictions using SARIMA(1,0,1)x(0,1,1)_12.
plot(m3.Temperatures,n1=as.vector(c(2000,1)),n.ahead=120,xlab='Year',type='o',
ylab='Temperatures',
main = "Predictions and prediction limits for the solar exposure
over the SARIMA(1,0,1)x(0,1,1)_12")
par(mfrow=c(1,1))
predict(m3.Temperatures,n.ahead = 24)## $pred
## Jan Feb Mar Apr May Jun Jul
## 2013
## 2014 19.728079 19.874177 17.805266 14.380814 11.537408 9.236016 8.632749
## 2015 19.709946 19.861503 17.796408 14.374623 11.533080 9.232992 8.630635
## Aug Sep Oct Nov Dec
## 2013 13.165444 15.974409 17.599915
## 2014 9.546103 11.225994 13.112340 15.937292 17.573972
## 2015 9.544625 11.224961
##
## $se
## Jan Feb Mar Apr May Jun Jul
## 2013
## 2014 1.0282584 1.0306934 1.0318809 1.0324606 1.0327436 1.0328818 1.0329494
## 2015 1.0354750 1.0354898 1.0354971 1.0355007 1.0355024 1.0355033 1.0355037
## Aug Sep Oct Nov Dec
## 2013 0.9914868 1.0129404 1.0232574
## 2014 1.0329824 1.0329985 1.0352560 1.0353837 1.0354461
## 2015 1.0355039 1.0355040As you can see, the seasonal component of the SARIMA has been captured by the pattern of Melbourne’s seasons. Further these predictions are quite accurate when looking at the years of 2001 to 2006. However, among the 14 years plotted from the data, most of the summers have higher average temperatures as compared to the future predictions made by the model. This could be a simple coincidence due to noise, or this could be due to a trend that this study failed to capture. Finally, this could be due to variability in the variance (heteroscedacity).
CONCLUSION
After careful analysis of Melbourne’s monthly temperatures between 1950 and 2013, the most accurate model for predictions was found to be the SARIMA(1,0,1)x(0,1,1)_12.
The predictions seem to be plausible on average, as the seasonal pattern is well captured by the SARIMA model. Nevertheless, one could argue that these predicted values do not exactly match the behaviour of Melbourne’s temperatures in preceeding years, particularly by having cooler summers.
In spite of excellent results on the overall model fit, it seems as if the time series shows a small change in behavior since the early nineties. We speculate that this is due to the appearance of a new trend, or to heteroscedasticity.