Forecasting Assignment 2
Student Details
- Shipren Jayadev (s3744421)
Load packages
library(readr)
library(x12)
library(dynlm)
library(dLagM)
library(forecast)
library(AER)
library(ggplot2)
library(TSA)
library(Hmisc)
library(tseries)
library(car)Load data
#Task 1
data1 <- read_csv("data1.csv")
data.x <- read_csv("data.x.csv")
#Task 2
data2 <- read_csv("data2.csv")Task 1
Introduction
In this task the goal is to determine which is the best model to be applied in a 2 year ahead forecast of solar radiation in a particular area based on monthly precipitation.
Data
The dataset being used consists of the monthly average horizontal solar radiation and the monthly precipitation series measured at the same points between January 1960 and December 2014.
Processing and visual desciption of data
solar <- as.vector(data1$solar)
ppt <- as.vector(data1$ppt)
solar.ts <- ts(solar, start = c(1960,1), end = c(2014,12), frequency = 12)
ppt.ts <- ts(ppt, start = c(1960,1), end = c(2014,12), frequency = 12)
plot(solar.ts)
plot(ppt.ts)
par(mfrow = c(1,2))
acf(solar.ts)
pacf(solar.ts)
acf(ppt.ts)
pacf(ppt.ts)
From the visualisations we can see the existence of seasonality in both figures with no obvious trends. There does appear to be an intervention point around the year 2000 in the solar time series as the J, A, M and S spikes are no longer as high. The same can be said for the PPT time series as around the year 2002 to 2004 there is a large spike in value.
From the acf and pacf of both variables we can also see that there is no apparent trend, however for both variables a seasonal pattern can be observed. As there is seasonality, the variables need to be adjusted to remove seasonality to determine if the time series is stationary
Transformations and Seasonality Adjustments
#Decomposition
solar.decom <- stl(solar.ts, t.window = 15, s.window = "periodic", robust = TRUE)
plot(solar.decom, main = "Solar Decomposition")
ppt.decom <- stl(ppt.ts, t.window = 15, s.window = "periodic", robust = TRUE)
plot(ppt.decom, main = "PPT Decomposition")
#Seasonally Adjusted Plots
solar.season.adjusted <- seasadj(solar.decom)
plot(solar.season.adjusted)
ppt.season.adjusted <- seasadj(ppt.decom)
plot(ppt.season.adjusted)
#Checking for stationarity
adf.test(solar.season.adjusted)## Warning in adf.test(solar.season.adjusted): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: solar.season.adjusted
## Dickey-Fuller = -6.1528, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationaryadf.test(ppt.season.adjusted)## Warning in adf.test(ppt.season.adjusted): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: ppt.season.adjusted
## Dickey-Fuller = -5.7359, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationaryppt.diff <- diff(ppt.season.adjusted)
plot(ppt.diff)
solar.final <- solar.season.adjusted
ppt.final <- ppt.season.adjustedFrom the adf tests, both time series have values less than 0.05 so we can conclude that they are both stationary after applying the differencing transformation.
The data is then decomposed and adjusted for seasonality producing the final variables we will be using for model fitting.
Lag Model Fitting
#Distributed lag model
print("MASE Values for dlm model:")## [1] "MASE Values for dlm model:"for ( i in 1:10){
model.dlm = dlm( x = as.vector(ppt.ts) , y = as.vector(solar.ts), q = i)
cat("q = ", i, "MASE = ", as.numeric(MASE(model.dlm)), "\n")
}## q = 1 MASE = 1.688457
## q = 2 MASE = 1.675967
## q = 3 MASE = 1.662703
## q = 4 MASE = 1.646357
## q = 5 MASE = 1.613848
## q = 6 MASE = 1.607532
## q = 7 MASE = 1.607042
## q = 8 MASE = 1.604806
## q = 9 MASE = 1.593121
## q = 10 MASE = 1.577996#When q = 10, MASE value is 1.578 (lowest)
model.dlm = dlm( x = as.vector(ppt.ts) , y = as.vector(solar.ts), q = 10)
summary(model.dlm)##
## Call:
## lm(formula = model.formula, data = design)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.9353 -5.4124 -0.7911 4.0184 30.8900
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.0105 1.0942 17.374 < 2e-16 ***
## x.t -7.3843 1.8995 -3.887 0.000112 ***
## x.1 -0.4763 2.5395 -0.188 0.851288
## x.2 -0.1324 2.5734 -0.051 0.958980
## x.3 1.7902 2.5781 0.694 0.487691
## x.4 1.9686 2.5808 0.763 0.445877
## x.5 3.4928 2.5807 1.353 0.176402
## x.6 0.5243 2.5787 0.203 0.838943
## x.7 1.6762 2.5797 0.650 0.516088
## x.8 0.9282 2.5673 0.362 0.717817
## x.9 0.3754 2.5338 0.148 0.882272
## x.10 -5.3798 1.8760 -2.868 0.004272 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.256 on 638 degrees of freedom
## Multiple R-squared: 0.3081, Adjusted R-squared: 0.2962
## F-statistic: 25.82 on 11 and 638 DF, p-value: < 2.2e-16
##
## AIC and BIC values for the model:
## AIC BIC
## 1 4602.658 4660.858checkresiduals(model.dlm$model$residuals)
In the chosen lag model, most of the coefficients are insignificant with the overall model being significant. The adjusted R-squared is low.
From the residuals check we can see that the histogram of residuals is symmetric and from the ACF we can see the residuals correlation is significant.
Dynamic Liner Model Fitting
Y.t <- solar.ts
T <- 480
S.t <- 1*(seq(Y.t) == T)
S.t.1 <- Lag(S.t,+1)
model.dynlm1 <- dynlm(Y.t ~ L(Y.t , k = 1 ) + S.t + trend(Y.t) + season(Y.t))
model.dynlm2 <- dynlm(Y.t ~ L(Y.t , k = 2 ) + S.t + trend(Y.t) + season(Y.t))
model.dynlm3 <- dynlm(Y.t ~ L(Y.t , k = 1 ) + S.t + S.t.1 + trend(Y.t) + season(Y.t))
model.dynlm4 <- dynlm(Y.t ~ L(Y.t , k = 1 ) + S.t + season(Y.t))
model.dynlm5 <- dynlm(Y.t ~ L(Y.t , k = 1 ) + S.t + trend(Y.t))
mase <- MASE(lm(model.dynlm1), lm(model.dynlm2), lm(model.dynlm3), lm(model.dynlm4), lm(model.dynlm5))
print(mase)## n MASE
## lm(model.dynlm1) 659 0.3740856
## lm(model.dynlm2) 658 0.5864478
## lm(model.dynlm3) 659 0.3737107
## lm(model.dynlm4) 659 0.3740093
## lm(model.dynlm5) 659 0.9729028#model.dynlm1 has the lowest MASE value of 0.375
summary(model.dynlm1)##
## Time series regression with "ts" data:
## Start = 1960(2), End = 2014(12)
##
## Call:
## dynlm(formula = Y.t ~ L(Y.t, k = 1) + S.t + trend(Y.t) + season(Y.t))
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8504 -0.9890 -0.1129 0.8326 16.8279
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2115296 0.3643660 3.325 0.000934 ***
## L(Y.t, k = 1) 0.9324176 0.0142368 65.494 < 2e-16 ***
## S.t 0.5993008 2.3191462 0.258 0.796170
## trend(Y.t) -0.0007248 0.0056472 -0.128 0.897917
## season(Y.t)Feb 1.9262264 0.4401138 4.377 1.41e-05 ***
## season(Y.t)Mar 4.8862999 0.4423924 11.045 < 2e-16 ***
## season(Y.t)Apr 3.7891154 0.4566473 8.298 6.23e-16 ***
## season(Y.t)May 5.1498542 0.4745335 10.852 < 2e-16 ***
## season(Y.t)Jun 3.3679048 0.5055095 6.662 5.79e-11 ***
## season(Y.t)Jul 1.6289566 0.5269001 3.092 0.002077 **
## season(Y.t)Aug -2.1607881 0.5344265 -4.043 5.91e-05 ***
## season(Y.t)Sep -4.7009039 0.5120670 -9.180 < 2e-16 ***
## season(Y.t)Oct -5.8126080 0.4781624 -12.156 < 2e-16 ***
## season(Y.t)Nov -5.0970411 0.4515972 -11.287 < 2e-16 ***
## season(Y.t)Dec -2.8434712 0.4432178 -6.416 2.72e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.297 on 644 degrees of freedom
## Multiple R-squared: 0.9463, Adjusted R-squared: 0.9452
## F-statistic: 811 on 14 and 644 DF, p-value: < 2.2e-16checkresiduals(model.dynlm1)
##
## Breusch-Godfrey test for serial correlation of order up to 24
##
## data: Residuals
## LM test = 132.76, df = 24, p-value < 2.2e-16From the selected model we can see that all coefficients except S.t and trend are significant.
Exponential Smoothing Fitting
#MASE is 0.6368203
model.ses <- ses(solar.ts)
summary(model.ses)##
## Forecast method: Simple exponential smoothing
##
## Model Information:
## Simple exponential smoothing
##
## Call:
## ses(y = solar.ts)
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 5.0575
##
## sigma: 4.576
##
## AIC AICc BIC
## 6296.371 6296.407 6309.847
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0001378357 4.569082 3.876391 -5.213129 27.30052 0.6368203
## ACF1
## Training set 0.6678374
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 5.148434 -0.7159726 11.01284 -3.820402 14.11727
## Feb 2015 5.148434 -3.1446744 13.44154 -7.534781 17.83165
## Mar 2015 5.148434 -5.0083386 15.30521 -10.385009 20.68188
## Apr 2015 5.148434 -6.5794989 16.87637 -12.787891 23.08476
## May 2015 5.148434 -7.9637280 18.26060 -14.904887 25.20175
## Jun 2015 5.148434 -9.2151719 19.51204 -16.818805 27.11567
## Jul 2015 5.148434 -10.3659966 20.66286 -18.578840 28.87571
## Aug 2015 5.148434 -11.4371603 21.73403 -20.217043 30.51391
## Sep 2015 5.148434 -12.4432208 22.74009 -21.755680 32.05255
## Oct 2015 5.148434 -13.3947777 23.69165 -23.210961 33.50783#MASE is 0.24716
model.hw1 <- hw(solar.ts,seasonal="additive")
summary(model.hw1)##
## Forecast method: Holt-Winters' additive method
##
## Model Information:
## Holt-Winters' additive method
##
## Call:
## hw(y = solar.ts, seasonal = "additive")
##
## Smoothing parameters:
## alpha = 0.9993
## beta = 0.0019
## gamma = 1e-04
##
## Initial states:
## l = 11.3306
## b = 0.209
## s = -10.7542 -8.4968 -3.227 3.0768 7.9264 10.9122
## 10.1422 7.3322 2.2353 -1.9528 -7.3626 -9.8316
##
## sigma: 2.3575
##
## AIC AICc BIC
## 5434.708 5435.662 5511.076
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.1179476 2.328736 1.504489 -2.230359 12.68101 0.24716 0.1703355
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 6.127399 3.1061581 9.148641 1.506810 10.74799
## Feb 2015 8.655602 4.3802222 12.930983 2.116973 15.19423
## Mar 2015 14.122221 8.8814747 19.362968 6.107191 22.13725
## Apr 2015 18.366625 12.3095941 24.423656 9.103196 27.63005
## May 2015 23.522264 16.7439512 30.300578 13.155729 33.88880
## Jun 2015 26.390991 18.9586817 33.823300 15.024255 37.75773
## Jul 2015 27.219228 19.1837595 35.254697 14.930039 39.50842
## Aug 2015 24.290519 15.6920179 32.889020 11.140246 37.44079
## Sep 2015 19.495874 10.3670365 28.624711 5.534522 33.45723
## Oct 2015 13.253765 3.6218809 22.885648 -1.476930 27.98446
## Nov 2015 8.042144 -2.0695730 18.153861 -7.422393 23.50668
## Dec 2015 5.840110 -4.7313921 16.411612 -10.327607 22.00783
## Jan 2016 6.819616 -4.1942247 17.833457 -10.024600 23.66383
## Feb 2016 9.347819 -2.0927719 20.788410 -8.149055 26.84469
## Mar 2016 14.814438 2.9609172 26.667959 -3.313958 32.94283
## Apr 2016 19.058842 6.8048110 31.312872 0.317919 37.79976
## May 2016 24.214481 11.5711780 36.857784 4.878218 43.55074
## Jun 2016 27.083208 14.0608586 40.105557 7.167243 46.99917
## Jul 2016 27.911445 14.5194061 41.303485 7.430089 48.39280
## Aug 2016 24.982736 11.2296053 38.735867 3.949138 46.01633
## Sep 2016 20.188091 6.0818048 34.294377 -1.385612 41.76179
## Oct 2016 13.945981 -0.5061083 28.398071 -8.156582 36.04855
## Nov 2016 8.734361 -6.0566988 23.525420 -13.886613 31.35533
## Dec 2016 6.532327 -8.5913309 21.655984 -16.597311 29.66196#MASE is 0.2461797
model.hw2 <- hw(solar.ts,seasonal="additive",damped = TRUE)
summary(model.hw2)##
## Forecast method: Damped Holt-Winters' additive method
##
## Model Information:
## Damped Holt-Winters' additive method
##
## Call:
## hw(y = solar.ts, seasonal = "additive", damped = TRUE)
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 1e-04
## gamma = 1e-04
## phi = 0.9388
##
## Initial states:
## l = 11.154
## b = 0.7632
## s = -10.4919 -8.137 -3.348 2.5794 8.08 11.1219
## 9.9586 6.9916 1.9573 -1.8565 -7.1607 -9.6946
##
## sigma: 2.3446
##
## AIC AICc BIC
## 5428.422 5429.489 5509.282
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01091357 2.314163 1.498521 -1.468083 12.44796 0.2461797
## ACF1
## Training set 0.1700724
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 5.945198 2.940530 8.949866 1.3499548 10.54044
## Feb 2015 8.479778 4.230548 12.729009 1.9811413 14.97842
## Mar 2015 13.784309 8.579902 18.988716 5.8248548 21.74376
## Apr 2015 17.598568 11.588776 23.608360 8.4073847 26.78975
## May 2015 22.632258 15.912803 29.351713 12.3557383 32.90878
## Jun 2015 25.599150 18.238024 32.960277 14.3412786 36.85702
## Jul 2015 26.761775 18.810497 34.713053 14.6013444 38.92221
## Aug 2015 23.722251 15.221610 32.222892 10.7216429 36.72286
## Sep 2015 18.222624 9.205955 27.239293 4.4328191 32.01243
## Oct 2015 12.293259 2.788471 21.798047 -2.2430606 26.82958
## Nov 2015 7.504219 -2.464877 17.473315 -7.7421977 22.75064
## Dec 2015 5.150488 -5.262286 15.563263 -10.7744758 21.07545
## Jan 2016 5.947275 -4.891175 16.785726 -10.6287044 22.52326
## Feb 2016 8.481728 -2.766251 19.729707 -8.7205712 25.68403
## Mar 2016 13.786139 2.142988 25.429291 -4.0205242 31.59280
## Apr 2016 17.600287 5.574906 29.625668 -0.7909464 35.99152
## May 2016 22.633871 10.238009 35.029734 3.6760353 41.59171
## Jun 2016 25.600665 12.845046 38.356283 6.0926298 45.10870
## Jul 2016 26.763197 13.657667 39.868726 6.7200186 46.80637
## Aug 2016 23.723586 10.277223 37.169950 3.1591478 44.28802
## Sep 2016 18.223877 4.445085 32.002669 -2.8489663 39.29672
## Oct 2016 12.294435 -1.808972 26.397843 -9.2748652 33.86374
## Nov 2016 7.505324 -6.915414 21.926061 -14.5492911 29.55994
## Dec 2016 5.151525 -9.579726 19.882776 -17.3779790 27.68103#MASE is 0.2233077
model.hw3 <- hw(solar.ts,seasonal="multiplicative")
summary(model.hw3)##
## Forecast method: Holt-Winters' multiplicative method
##
## Model Information:
## Holt-Winters' multiplicative method
##
## Call:
## hw(y = solar.ts, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha = 0.8251
## beta = 1e-04
## gamma = 0.0259
##
## Initial states:
## l = 10.5659
## b = -0.2371
## s = 0.5024 0.6292 0.9319 1.2229 1.4923 1.6309
## 1.5606 1.3672 1.0278 0.8128 0.5106 0.3114
##
## sigma: 0.3971
##
## AIC AICc BIC
## 6648.746 6649.699 6725.114
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.04281136 2.12539 1.359297 -0.4896895 10.87401 0.2233077
## ACF1
## Training set 0.06551473
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 5.610335 2.75534413 8.465326 1.2440033 9.976666
## Feb 2015 6.512806 2.03809262 10.987520 -0.3306776 13.356290
## Mar 2015 8.786120 1.33667747 16.235562 -2.6068190 20.179058
## Apr 2015 10.192785 -0.02672876 20.412298 -5.4366123 25.822181
## May 2015 12.512597 -1.96157245 26.986766 -9.6237347 34.648928
## Jun 2015 14.061609 -4.41048190 32.533701 -14.1890164 42.312235
## Jul 2015 14.502088 -6.89909775 35.903274 -18.2282011 47.232377
## Aug 2015 12.945620 -8.35024040 34.241481 -19.6235881 45.514829
## Sep 2015 10.352210 -8.52365174 29.228072 -18.5159294 39.220350
## Oct 2015 7.418279 -7.51119146 22.347750 -15.4143760 30.250935
## Nov 2015 4.989730 -6.05900058 16.038460 -11.9078451 21.887305
## Dec 2015 3.871776 -5.53876625 13.282319 -10.5204066 18.263960
## Jan 2016 4.199400 -7.01753622 15.416336 -12.9554236 21.354224
## Feb 2016 4.838892 -9.30305408 18.980837 -16.7893479 26.467131
## Mar 2016 6.477174 -14.21845674 27.172805 -25.1740619 38.128410
## Apr 2016 7.452633 -18.56883703 33.474103 -32.3437711 47.249037
## May 2016 9.069749 -25.52978297 43.669280 -43.8456686 61.985166
## Jun 2016 10.099485 -31.99838423 52.197355 -54.2836502 74.482621
## Jul 2016 10.315199 -36.68017471 57.310572 -61.5580226 82.188420
## Aug 2016 9.113765 -36.29212001 54.519651 -60.3285438 78.556075
## Sep 2016 7.208700 -32.09306283 46.510463 -52.8981594 67.315559
## Oct 2016 5.105869 -25.38360449 35.595343 -41.5237568 51.735495
## Nov 2016 3.391945 -18.81647619 25.600367 -30.5729044 37.356795
## Dec 2016 2.597256 -16.07148862 21.266001 -25.9541251 31.148637#MASE is 0.2320404
model.hw4 <- hw(solar.ts,seasonal="multiplicative",exponential = TRUE)
summary(model.hw4)##
## Forecast method: Holt-Winters' multiplicative method with exponential trend
##
## Model Information:
## Holt-Winters' multiplicative method with exponential trend
##
## Call:
## hw(y = solar.ts, seasonal = "multiplicative", exponential = TRUE)
##
## Smoothing parameters:
## alpha = 0.6499
## beta = 1e-04
## gamma = 0.0597
##
## Initial states:
## l = 9.3759
## b = 0.9889
## s = 0.3552 0.4789 0.952 1.0553 1.4608 1.7926
## 1.6746 1.4529 1.1145 0.8672 0.5333 0.2627
##
## sigma: 0.3755
##
## AIC AICc BIC
## 6584.208 6585.161 6660.576
##
## Error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.06388152 2.208727 1.412454 -1.303792 11.28449 0.2320404 0.153871
##
## Forecasts:
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 5.797870 2.9001990 8.530608 1.4149978 9.88812
## Feb 2015 7.513359 3.2834470 12.021152 1.6873557 14.91881
## Mar 2015 10.758212 4.3148647 18.322565 2.1012864 24.19686
## Apr 2015 12.694087 4.5843742 22.648720 2.3635223 31.71288
## May 2015 15.340074 4.8910884 28.292417 2.4597857 40.50073
## Jun 2015 17.239431 4.8327553 32.760516 2.1893904 49.80672
## Jul 2015 18.004700 4.7053905 35.022252 2.2886088 54.10633
## Aug 2015 16.206596 4.0144303 32.386812 1.9099695 49.65613
## Sep 2015 13.011256 2.9221156 26.272001 1.2919208 41.76711
## Oct 2015 9.215869 1.9066318 18.907877 0.8356363 30.08522
## Nov 2015 6.214762 1.1901737 12.955522 0.4866940 22.17772
## Dec 2015 4.565667 0.8013365 9.810552 0.3309239 17.44177
## Jan 2016 5.220661 0.8074575 11.256801 0.3273246 20.84995
## Feb 2016 6.765364 1.0143805 15.253213 0.3742710 28.39995
## Mar 2016 9.687175 1.3329042 21.925952 0.5147698 41.27373
## Apr 2016 11.430324 1.5014058 25.790960 0.5811961 54.95675
## May 2016 13.812889 1.6551989 31.689776 0.6390718 63.31986
## Jun 2016 15.523155 1.6963109 35.517387 0.6312494 74.63215
## Jul 2016 16.212237 1.6735103 38.322188 0.5959797 80.57812
## Aug 2016 14.593143 1.3856435 34.628247 0.4955130 73.36872
## Sep 2016 11.715916 1.0527904 27.360447 0.3795664 60.77175
## Oct 2016 8.298381 0.6880969 19.919345 0.2501803 44.34245
## Nov 2016 5.596049 0.4380251 13.156356 0.1700061 29.07893
## Dec 2016 4.111130 0.3112401 9.669491 0.1001913 21.60976#Lowest MASE value belongs to model.hw3
checkresiduals(model.hw3)
##
## Ljung-Box test
##
## data: Residuals from Holt-Winters' multiplicative method
## Q* = 38.585, df = 8, p-value = 5.868e-06
##
## Model df: 16. Total lags used: 24From the model residuals we can see that there are 2 significant points in the ACF and the histogram of residuals is mostly symmetric with a few outliers. However we can still see some seasonality in the ACF plot.
State Space Model Fitting
#MASE is 0.4334691
model1.ets <- ets(solar.ts, model = "AAN")
summary(model1.ets)## ETS(A,Ad,N)
##
## Call:
## ets(y = solar.ts, model = "AAN")
##
## Smoothing parameters:
## alpha = 0.9279
## beta = 0.9279
## phi = 0.8
##
## Initial states:
## l = 13.7196
## b = 2.5786
##
## sigma: 3.4657
##
## AIC AICc BIC
## 5932.524 5932.653 5959.478
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.008312543 3.452593 2.638571 3.790894 19.25631 0.4334691
## ACF1
## Training set 0.07272876#MASE is 0.254704
model2.ets <- ets(solar.ts, model = "ANA")
summary(model2.ets)## ETS(A,N,A)
##
## Call:
## ets(y = solar.ts, model = "ANA")
##
## Smoothing parameters:
## alpha = 0.9997
## gamma = 3e-04
##
## Initial states:
## l = 22.6142
## s = -10.2887 -8.5269 -4.0534 1.7184 7.355 10.8208
## 10.4774 7.7101 2.4749 -1.7156 -6.7936 -9.1783
##
## sigma: 2.3884
##
## AIC AICc BIC
## 5449.974 5450.719 5517.358
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0108645 2.362974 1.55041 -2.640142 12.95336 0.254704 0.1825724#MASE is 0.6513597
model3.ets <- ets(solar.ts, model = "MAN")
summary(model3.ets)## ETS(M,A,N)
##
## Call:
## ets(y = solar.ts, model = "MAN")
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.0018
##
## Initial states:
## l = 7.8559
## b = 2.4488
##
## sigma: 0.3004
##
## AIC AICc BIC
## 6433.364 6433.456 6455.825
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.466339 4.822041 3.964894 -17.35094 30.69019 0.6513597 0.6697971#MASE is 0.2461797
model4.ets <- ets(solar.ts, model = "AAA")
summary(model4.ets)## ETS(A,Ad,A)
##
## Call:
## ets(y = solar.ts, model = "AAA")
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 1e-04
## gamma = 1e-04
## phi = 0.9388
##
## Initial states:
## l = 11.154
## b = 0.7632
## s = -10.4919 -8.137 -3.348 2.5794 8.08 11.1219
## 9.9586 6.9916 1.9573 -1.8565 -7.1607 -9.6946
##
## sigma: 2.3446
##
## AIC AICc BIC
## 5428.422 5429.489 5509.282
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01091357 2.314163 1.498521 -1.468083 12.44796 0.2461797
## ACF1
## Training set 0.1700724checkresiduals(model4.ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 210.76, df = 7, p-value < 2.2e-16
##
## Model df: 17. Total lags used: 24#MASE is 0.695816
model5.ets <- ets(solar.ts, model = "MMN")
summary(model5.ets)## ETS(M,M,N)
##
## Call:
## ets(y = solar.ts, model = "MMN")
##
## Smoothing parameters:
## alpha = 0.9999
## beta = 0.0063
##
## Initial states:
## l = 9.1784
## b = 1.1871
##
## sigma: 0.3511
##
## AIC AICc BIC
## 6609.163 6609.255 6631.624
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.644301 5.221518 4.235504 -14.9362 30.65972 0.695816 0.6833568#MASE is 0.3222574
model6.ets <- ets(solar.ts, model = "MAM")
summary(model6.ets)## ETS(M,Ad,M)
##
## Call:
## ets(y = solar.ts, model = "MAM")
##
## Smoothing parameters:
## alpha = 0.7922
## beta = 6e-04
## gamma = 0.0707
## phi = 0.98
##
## Initial states:
## l = 9.3693
## b = 1.2404
## s = 0.6953 0.3066 0.5802 0.9879 1.3594 1.5469
## 1.4723 1.4055 1.227 0.9573 0.6959 0.7656
##
## sigma: 0.2274
##
## AIC AICc BIC
## 5953.502 5954.569 6034.363
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.1750605 3.02791 1.961614 -5.497209 17.21119 0.3222574 0.2565219#MASE is 0.6368203
model7.ets <- ets(solar.ts, model = "ANN")
summary(model7.ets)## ETS(A,N,N)
##
## Call:
## ets(y = solar.ts, model = "ANN")
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 5.0575
##
## sigma: 4.576
##
## AIC AICc BIC
## 6296.371 6296.407 6309.847
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0001378357 4.569082 3.876391 -5.213129 27.30052 0.6368203
## ACF1
## Training set 0.6678374#MASE is 0.6369599
model8.ets <- ets(solar.ts, model = "MNN")
summary(model8.ets)## ETS(M,N,N)
##
## Call:
## ets(y = solar.ts, model = "MNN")
##
## Smoothing parameters:
## alpha = 0.9999
##
## Initial states:
## l = 4.4856
##
## sigma: 0.3868
##
## AIC AICc BIC
## 6619.776 6619.812 6633.253
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.001004352 4.569136 3.877241 -5.195978 27.31733 0.6369599
## ACF1
## Training set 0.6678785The model with the lowest MASE value is model4.ets with model set to AAA.
From the residuals we can see there are a few significant values in the ACF. And the histogram of residuals is symmetric.
Forecasting
From the above model fitting, the model with the lowest MASE value is model.hw3 which uses the holt-winters model with seasonal set to multiplicative. model.hw3 will be used to predict the 2 year ahead forecast.
forecast(model.hw3)## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Jan 2015 5.610335 2.75534413 8.465326 1.2440033 9.976666
## Feb 2015 6.512806 2.03809262 10.987520 -0.3306776 13.356290
## Mar 2015 8.786120 1.33667747 16.235562 -2.6068190 20.179058
## Apr 2015 10.192785 -0.02672876 20.412298 -5.4366123 25.822181
## May 2015 12.512597 -1.96157245 26.986766 -9.6237347 34.648928
## Jun 2015 14.061609 -4.41048190 32.533701 -14.1890164 42.312235
## Jul 2015 14.502088 -6.89909775 35.903274 -18.2282011 47.232377
## Aug 2015 12.945620 -8.35024040 34.241481 -19.6235881 45.514829
## Sep 2015 10.352210 -8.52365174 29.228072 -18.5159294 39.220350
## Oct 2015 7.418279 -7.51119146 22.347750 -15.4143760 30.250935
## Nov 2015 4.989730 -6.05900058 16.038460 -11.9078451 21.887305
## Dec 2015 3.871776 -5.53876625 13.282319 -10.5204066 18.263960
## Jan 2016 4.199400 -7.01753622 15.416336 -12.9554236 21.354224
## Feb 2016 4.838892 -9.30305408 18.980837 -16.7893479 26.467131
## Mar 2016 6.477174 -14.21845674 27.172805 -25.1740619 38.128410
## Apr 2016 7.452633 -18.56883703 33.474103 -32.3437711 47.249037
## May 2016 9.069749 -25.52978297 43.669280 -43.8456686 61.985166
## Jun 2016 10.099485 -31.99838423 52.197355 -54.2836502 74.482621
## Jul 2016 10.315199 -36.68017471 57.310572 -61.5580226 82.188420
## Aug 2016 9.113765 -36.29212001 54.519651 -60.3285438 78.556075
## Sep 2016 7.208700 -32.09306283 46.510463 -52.8981594 67.315559
## Oct 2016 5.105869 -25.38360449 35.595343 -41.5237568 51.735495
## Nov 2016 3.391945 -18.81647619 25.600367 -30.5729044 37.356795
## Dec 2016 2.597256 -16.07148862 21.266001 -25.9541251 31.148637plot(forecast(model.hw3), main = "2 Year Ahead Forecast of Solar Radiation", xlab = "Year", ylab = "Radiation")
TASK 2
Introduction
The objective of this task is to determine if the correlation between residential property price index and quarterly population changes is spurious or not.
The Data
The data contains the quarterly property price index(PPI) and changes in the population in Melbourne between September 2003 and December 2016.
Converting to time series and visualisation
#Converting to time series
ppi.ts <- ts(data2[, 2:3], start = c(2003), frequency = 4)
plot(ppi.ts)
From the first plot we can see that both patterns are increasing in roughly the same way, implying correlation.
Checking CCF
ccf(as.vector(ppi.ts[,1]), as.vector(ppi.ts[,2]))
The majority of cross-correlations are significant so we will need to apply prewhitening to determine if it is spurious or not.
Prewhitening
me.dif <- ts.intersect(diff(diff(ppi.ts[,1],4)),
diff(diff(ppi.ts[,2],4)))
prewhiten(as.vector(me.dif[,1]),as.vector(me.dif[,2]),ylab='CCF', main="CFF after prewhitening ")
Looking at the CCF now we can see that the correlation is spurious after applying the prewhitening process.