MATH1307 Forecasting- Assignment 2
Syed Wajahath - s3750039
24/09/2019
R Markdown
rm(list=ls()) library(dynlm) library(ggplot2) library(AER) library(Hmisc) library(forecast) library(x12) library(dLagM) library(TSA) library(readr) library(plyr) library(MuMIn)
Task 1
AIM: The aim of this task is to forecast the quantam of horizontal solar radiation falling on ground at a particular location across the world. To accomplish the task, related data of solar radiation and precipitation aggregated monthly over from January 1960 to December 2014 is used. Forecasting here is done using related precipitation data based on MASE using dLagM package, dynlm package, exponential smoothing and corresponding state-space models.
Data
Loading the given data to make it a series and further auto-scaling it to plot.
setwd("C:/Users/Wajahath/Desktop/Analytics/Sem 2/Forcasting/Assignmet 2")
solradppt <- read.csv("data1.csv")
solrad <- ts(solradppt$solar, start = c(1960,1), frequency = 12)
ppt <- ts(solradppt$ppt, start= c(1960,1), frequency = 12)
solppt <- ts.intersect(solrad, ppt) #Joining time series of radiation and precipitation
forecastdata <- read.csv("data.x.csv")
forecastdata_ts <- ts(forecastdata)
#Scaling the data as they are on different scales
scaled_solppt = scale(solppt)
par(oma=c(0,0,2,0))
plot(scaled_solppt, plot.type="s",col = c("steelblue4", "sienna1"),
main = "Solar Radiation vs Precipitation over years")
legend("topleft",lty=1, text.width = 11, col=c("steelblue4", "sienna1"), c("Solar Radiation", "Ppt"))
Descriptive Plots
- ACF and PACF plots
#ACF and PACF
par(mar=c(4,4.5,5,1),mfrow=c(1,2),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)
acf(solrad, main = "ACF plot for Solar Radiation series")
pacf(solrad, main = "PACF plot for Solar Radiation series")
From the ACF and PACF plots above,
The wavy pattern of ACF suggest the existence of seasonality in the series.
There is no evidence of changing variance.
Patterns of moving average is seen by existence of successive points.
The lengthier first line in PACF suggests the possibility of trend.
There seems to be no intervention at any point.
Now, let us cross check for the existence of seasonality using the decomposition method.
#Seasonality check
plotSeasFac(x12(solrad))
It is clear from the above plot how the solar radiation is lower in the months of January and December. The graphs depicts a curve which goes up and comes down in nearly uniform way.
Model Fiiting
To find which is the best suited model, a systematic approach needs to be followed. The foremost step is to check the corelation between the series. Following which we can fit the best model out of existing Distributed LM / Dynamic LM / Exponential smoothing / State-space models.
#Corelation test
cor(solrad,ppt)## [1] -0.4540277Distributed Lag Models
In order to have best results for the forecast, let us find which model fits the best. The different models used are;
- DLM
- Poly DLM
- Koyck DLM
- AR DLM
- DLM - Using DLM, trying to find the best suited significant model with initial q to be 1 till q = 10.
model_dlm <- data.frame(Model=character() , MASE=numeric() ,
BIC= numeric() , AICC=numeric() , AIC=numeric())
for ( i in 1:10){
dlm = dlm( y = as.vector(solrad) , x = as.vector(ppt), q = i )
model_dlm = rbind.fill(model_dlm,
cbind(Model= paste0("DLM q= ",i), MASE=MASE(dlm$model),
BIC=BIC(dlm$model),AICC = AICc(dlm$model),
AIC = AIC(dlm$model)))
}
model_dlm## Model MASE BIC AICC AIC
## 1 DLM q= 1 1.688457 4746.676 4728.774 4728.713
## 2 DLM q= 2 1.675967 4735.095 4712.741 4712.649
## 3 DLM q= 3 1.662703 4715.478 4688.681 4688.551
## 4 DLM q= 4 1.646357 4695.003 4663.772 4663.600
## 5 DLM q= 5 1.613848 4680.499 4644.845 4644.622
## 6 DLM q= 6 1.607532 4677.837 4637.769 4637.489
## 7 DLM q= 7 1.607042 4677.532 4633.059 4632.716
## 8 DLM q= 8 1.604806 4675.267 4626.399 4625.986
## 9 DLM q= 9 1.593121 4668.827 4615.573 4615.084
## 10 DLM q= 10 1.577996 4660.858 4603.230 4602.658- The best DLM interms of MASE, AIc, AICc and BIC is found to be when q = 10 with MASE value 1.5779, AIC value 4602.658, AICc value 4603.23 and BIC value 4660.858
To check other parameters and residuals, do Breusch-Godfrey test.
bgtest(dlm$model)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: dlm$model
## LM test = 522.18, df = 1, p-value < 2.2e-16checkresiduals(dlm$model$residuals, test = F)
- The BG tests shows significant auto-correation in data which deems this model to be not good. Also, the residual plot shows residuals not being distributed normally.
This can be confirmed by vif and Shapiro test.
vif(dlm$model)## x.t x.1 x.2 x.3 x.4 x.5 x.6 x.7
## 4.244594 7.665259 7.910115 7.952633 7.957841 7.941836 7.911213 7.901898
## x.8 x.9 x.10
## 7.847965 7.653512 4.228221- The results of vif and Shapiro tests suggest the DLM with q = 10 is not the best suited model which promts us to check with Poly DLm.
- Polynomial DLM - Using Ploy DLM, trying to find the best model with i = 1 to 10 and k =2.
model_dlm <- data.frame(Model=character(),MASE=numeric(),
BIC= numeric(),AICC=numeric(),AIC=numeric())
for ( i in 1:10){
polydlm = polyDlm( y = as.vector(solrad) , x = as.vector(ppt), q = i, k= 2, show.beta = FALSE )
model_dlm = rbind.fill(model_dlm,
cbind(Model= paste0("Poly DLM q = ",i, " k = 2" ), MASE=MASE(polydlm$model),
BIC=BIC(polydlm$model),
AICC =AICc(polydlm$model),
AIC = AIC(polydlm$model)))
}
model_dlm## Model MASE BIC AICC AIC
## 1 Poly DLM q = 1 k = 2 1.688457 4746.676 4728.774 4728.713
## 2 Poly DLM q = 2 k = 2 1.675967 4735.095 4712.741 4712.649
## 3 Poly DLM q = 3 k = 2 1.666349 4711.457 4689.111 4689.018
## 4 Poly DLM q = 4 k = 2 1.653286 4687.171 4664.833 4664.741
## 5 Poly DLM q = 5 k = 2 1.631939 4667.673 4645.342 4645.250
## 6 Poly DLM q = 6 k = 2 1.619987 4656.942 4634.619 4634.526
## 7 Poly DLM q = 7 k = 2 1.619091 4650.382 4628.067 4627.974
## 8 Poly DLM q = 8 k = 2 1.619034 4642.870 4620.563 4620.470
## 9 Poly DLM q = 9 k = 2 1.608805 4630.253 4607.954 4607.861
## 10 Poly DLM q = 10 k = 2 1.592555 4614.289 4591.997 4591.904- Koyck DLM
koyck = koyckDlm(y=as.vector(solrad),x=as.vector(ppt))
model_dlm <- data.frame(Model=character(),MASE=numeric(),
BIC= numeric(),AICC=numeric(),AIC=numeric())
attr(koyck$model, "class")="lm"
Koyckaic <- AIC(koyck$model)
model_dlm = rbind.fill(model_dlm,cbind(Model="Koyck",MASE=MASE(koyck),
BIC = NA,
AICC = NA,
AIC = Koyckaic ))
model_dlm## Model MASE BIC AICC AIC
## 1 Koyck 1.032483 NA NA 3946.4764.ARDLM
model_dlm <- data.frame(Model=character(),MASE=numeric(),
BIC= numeric(),AICC=numeric(),AIC=numeric())
for ( i in 1:10){
ardlm = ardlDlm( y = as.vector(solrad) , x = as.vector(ppt), p = i, q= i)
model_dlm = rbind.fill(model_dlm,
cbind(Model= paste0("ARDLM q = ",i, " k = ",i ), MASE=MASE(ardlm),
BIC=BIC(ardlm$model),
AICC =AICc(ardlm$model),
AIC = AIC(ardlm$model)))
}
model_dlm## Model MASE BIC AICC AIC
## 1 ARDLM q = 1 k = 1 0.8392434 3734.765 3712.403 3712.311
## 2 ARDLM q = 2 k = 2 0.4951319 3260.476 3229.224 3229.051
## 3 ARDLM q = 3 k = 3 0.4737144 3179.798 3139.687 3139.409
## 4 ARDLM q = 4 k = 4 0.4665123 3180.772 3131.834 3131.424
## 5 ARDLM q = 5 k = 5 0.4479311 3156.177 3098.445 3097.877
## 6 ARDLM q = 6 k = 6 0.4205664 3098.223 3031.728 3030.976
## 7 ARDLM q = 7 k = 7 0.4103932 3078.600 3003.377 3002.413
## 8 ARDLM q = 8 k = 8 0.4064415 3078.922 2995.003 2993.801
## 9 ARDLM q = 9 k = 9 0.4042283 3069.631 2977.052 2975.583
## 10 ARDLM q = 10 k = 10 0.3996089 3065.571 2964.364 2962.601ardl10 = ardlDlm(y=as.vector(solrad),x=as.vector(ppt),p = 10,q = 10)
checkresiduals(ardl10$model, test = F)
Dynamic Lag Models
Y.t = solrad
X.t = ppt
dynlm1 = dynlm(Y.t ~ L(Y.t ,1 ) + trend(Y.t) + season(Y.t))
summary(dynlm1)##
## Time series regression with "ts" data:
## Start = 1960(2), End = 2014(12)
##
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + trend(Y.t) + season(Y.t))
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8501 -0.9906 -0.1119 0.8325 16.8276
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.2105921 0.3640843 3.325 0.000934 ***
## L(Y.t, 1) 0.9323820 0.0142258 65.542 < 2e-16 ***
## trend(Y.t) -0.0006813 0.0056406 -0.121 0.903905
## season(Y.t)Feb 1.9262673 0.4397952 4.380 1.39e-05 ***
## season(Y.t)Mar 4.8864291 0.4420720 11.053 < 2e-16 ***
## season(Y.t)Apr 3.7894322 0.4563152 8.304 5.90e-16 ***
## season(Y.t)May 5.1503064 0.4741868 10.861 < 2e-16 ***
## season(Y.t)Jun 3.3685316 0.5051379 6.669 5.56e-11 ***
## season(Y.t)Jul 1.6296826 0.5265113 3.095 0.002052 **
## season(Y.t)Aug -2.1600321 0.5340317 -4.045 5.87e-05 ***
## season(Y.t)Sep -4.7002552 0.5116903 -9.186 < 2e-16 ***
## season(Y.t)Oct -5.8121501 0.4778130 -12.164 < 2e-16 ***
## season(Y.t)Nov -5.0968011 0.4512694 -11.294 < 2e-16 ***
## season(Y.t)Dec -2.8325126 0.4408651 -6.425 2.56e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.295 on 645 degrees of freedom
## Multiple R-squared: 0.9463, Adjusted R-squared: 0.9452
## F-statistic: 874.6 on 13 and 645 DF, p-value: < 2.2e-16dynlm2 = dynlm(Y.t ~ L(Y.t ,1 ) + L(Y.t ,2 ) + trend(Y.t) + season(Y.t))
summary(dynlm2)##
## Time series regression with "ts" data:
## Start = 1960(3), End = 2014(12)
##
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + trend(Y.t) + season(Y.t))
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.2013 -1.0035 -0.0807 0.8270 16.9555
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8122916 0.3719451 4.872 1.39e-06 ***
## L(Y.t, 1) 1.1319003 0.0385108 29.392 < 2e-16 ***
## L(Y.t, 2) -0.2141527 0.0385114 -5.561 3.94e-08 ***
## trend(Y.t) -0.0006783 0.0055291 -0.123 0.902394
## season(Y.t)Feb 1.3206473 0.4468083 2.956 0.003233 **
## season(Y.t)Mar 3.8878770 0.4681882 8.304 5.95e-16 ***
## season(Y.t)Apr 2.2727488 0.5231550 4.344 1.62e-05 ***
## season(Y.t)May 4.0033234 0.5077132 7.885 1.35e-14 ***
## season(Y.t)Jun 2.0598056 0.5474589 3.762 0.000184 ***
## season(Y.t)Jul 0.8170809 0.5354520 1.526 0.127510
## season(Y.t)Aug -2.5446968 0.5269572 -4.829 1.72e-06 ***
## season(Y.t)Sep -4.3021817 0.5054747 -8.511 < 2e-16 ***
## season(Y.t)Oct -4.9890230 0.4900800 -10.180 < 2e-16 ***
## season(Y.t)Nov -4.1995146 0.4698438 -8.938 < 2e-16 ***
## season(Y.t)Dec -2.2469226 0.4438227 -5.063 5.41e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.245 on 643 degrees of freedom
## Multiple R-squared: 0.9487, Adjusted R-squared: 0.9476
## F-statistic: 849.5 on 14 and 643 DF, p-value: < 2.2e-16dynlm3 = dynlm(Y.t ~ L(Y.t ,1 ) + L(Y.t ,2 ) + X.t + trend(Y.t) + season(Y.t))
summary(dynlm3)##
## Time series regression with "ts" data:
## Start = 1960(3), End = 2014(12)
##
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + X.t + trend(Y.t) +
## season(Y.t))
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.1317 -0.9985 -0.0884 0.8337 17.0093
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.983337 0.436286 4.546 6.54e-06 ***
## L(Y.t, 1) 1.131707 0.038525 29.376 < 2e-16 ***
## L(Y.t, 2) -0.214438 0.038526 -5.566 3.83e-08 ***
## X.t -0.243985 0.324990 -0.751 0.453081
## trend(Y.t) -0.000493 0.005536 -0.089 0.929069
## season(Y.t)Feb 1.299241 0.447869 2.901 0.003848 **
## season(Y.t)Mar 3.862038 0.469610 8.224 1.09e-15 ***
## season(Y.t)Apr 2.245544 0.524586 4.281 2.15e-05 ***
## season(Y.t)May 3.952718 0.512339 7.715 4.64e-14 ***
## season(Y.t)Jun 1.958115 0.564147 3.471 0.000553 ***
## season(Y.t)Jul 0.685603 0.563537 1.217 0.224201
## season(Y.t)Aug -2.680867 0.557469 -4.809 1.89e-06 ***
## season(Y.t)Sep -4.420256 0.529541 -8.347 4.29e-16 ***
## season(Y.t)Oct -5.038937 0.494734 -10.185 < 2e-16 ***
## season(Y.t)Nov -4.206824 0.470104 -8.949 < 2e-16 ***
## season(Y.t)Dec -2.221271 0.445286 -4.988 7.85e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.245 on 642 degrees of freedom
## Multiple R-squared: 0.9488, Adjusted R-squared: 0.9476
## F-statistic: 792.3 on 15 and 642 DF, p-value: < 2.2e-16dynlm4 = dynlm(Y.t ~ L(Y.t ,1 ) + L(Y.t ,2 ) +
X.t + L(X.t,1) + trend(Y.t) + season(Y.t))
summary(dynlm4)##
## Time series regression with "ts" data:
## Start = 1960(3), End = 2014(12)
##
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + L(Y.t, 2) + X.t + L(X.t, 1) +
## trend(Y.t) + season(Y.t))
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.1184 -1.0169 -0.1062 0.8286 16.8855
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.0433705 0.4559411 4.482 8.77e-06 ***
## L(Y.t, 1) 1.1310883 0.0385724 29.324 < 2e-16 ***
## L(Y.t, 2) -0.2139682 0.0385638 -5.548 4.22e-08 ***
## X.t -0.0692011 0.5023997 -0.138 0.890488
## L(X.t, 1) -0.2291443 0.5020643 -0.456 0.648254
## trend(Y.t) -0.0004628 0.0055403 -0.084 0.933450
## season(Y.t)Feb 1.2931622 0.4483428 2.884 0.004054 **
## season(Y.t)Mar 3.8426013 0.4718258 8.144 2.00e-15 ***
## season(Y.t)Apr 2.2243271 0.5269639 4.221 2.78e-05 ***
## season(Y.t)May 3.9468506 0.5128165 7.696 5.31e-14 ***
## season(Y.t)Jun 1.9674717 0.5648678 3.483 0.000529 ***
## season(Y.t)Jul 0.6672213 0.5653213 1.180 0.238338
## season(Y.t)Aug -2.7257717 0.5664231 -4.812 1.86e-06 ***
## season(Y.t)Sep -4.4859263 0.5490569 -8.170 1.64e-15 ***
## season(Y.t)Oct -5.1393686 0.5417433 -9.487 < 2e-16 ***
## season(Y.t)Nov -4.2751836 0.4936640 -8.660 < 2e-16 ***
## season(Y.t)Dec -2.2728364 0.4596625 -4.945 9.76e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.247 on 641 degrees of freedom
## Multiple R-squared: 0.9488, Adjusted R-squared: 0.9475
## F-statistic: 741.9 on 16 and 641 DF, p-value: < 2.2e-16dynlm5 = dynlm(Y.t ~ L(Y.t ,1 ) + season(Y.t))
summary(dynlm5)##
## Time series regression with "ts" data:
## Start = 1960(2), End = 2014(12)
##
## Call:
## dynlm(formula = Y.t ~ L(Y.t, 1) + season(Y.t))
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8554 -1.0002 -0.1205 0.8259 16.8225
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.19188 0.32924 3.620 0.000317 ***
## L(Y.t, 1) 0.93237 0.01421 65.592 < 2e-16 ***
## season(Y.t)Feb 1.92656 0.43945 4.384 1.36e-05 ***
## season(Y.t)Mar 4.88669 0.44173 11.063 < 2e-16 ***
## season(Y.t)Apr 3.78970 0.45596 8.311 5.58e-16 ***
## season(Y.t)May 5.15056 0.47382 10.870 < 2e-16 ***
## season(Y.t)Jun 3.36878 0.50475 6.674 5.36e-11 ***
## season(Y.t)Jul 1.62991 0.52611 3.098 0.002032 **
## season(Y.t)Aug -2.15985 0.53362 -4.048 5.80e-05 ***
## season(Y.t)Sep -4.70016 0.51130 -9.193 < 2e-16 ***
## season(Y.t)Oct -5.81217 0.47745 -12.173 < 2e-16 ***
## season(Y.t)Nov -5.09695 0.45092 -11.303 < 2e-16 ***
## season(Y.t)Dec -2.83277 0.44052 -6.430 2.47e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.293 on 646 degrees of freedom
## Multiple R-squared: 0.9463, Adjusted R-squared: 0.9453
## F-statistic: 948.9 on 12 and 646 DF, p-value: < 2.2e-16model_dlm <- data.frame(Model=character(),MASE=numeric(),
BIC= numeric(),AICC=numeric(),AIC=numeric())
model_dlm = rbind(model_dlm,cbind(Model="DYNLM 5",
BIC = BIC(dynlm5),
AICC = AICc(dynlm5),
AIC = AIC(dynlm5)))
model_dlm## Model BIC AICC AIC
## 1 DYNLM 5 3041.90343987038 2979.68548430039 2979.03331038734checkresiduals(dynlm5 ,test = F)
Exponential Smoothing Methods
model_dlm <- data.frame(Model=character() , MASE=numeric() ,
BIC= numeric() , AICC=numeric() , AIC=numeric())
hw_addi = hw(solrad,seasonal="additive")
hw_multip = hw(solrad,seasonal="multiplicative")
hw_multipDmp = hw(solrad,seasonal="multiplicative" , damped = T)
hw_multiDmpExp = hw(solrad,seasonal="multiplicative" , damped = T, exponential = T)
hw_multipExp = hw(solrad,seasonal="multiplicative" , damped = F, exponential = T)
model_dlm = rbind(model_dlm,cbind(Model="hw_addi",MASE=accuracy(hw_addi)[6],
BIC = hw_addi$model$bic,
AICC =hw_addi$model$aicc,
AIC = hw_addi$model$aic))
model_dlm = rbind(model_dlm,cbind(Model="hw_multip",MASE=accuracy(hw_multip)[6],
BIC = hw_multip$model$bic,
AICC = hw_multip$model$aicc,
AIC = hw_multip$model$aic))
model_dlm = rbind(model_dlm,cbind(Model="hw_multipDmp",MASE=accuracy(hw_multipDmp)[6],
BIC = hw_multipDmp$model$bic,
AICC = hw_multipDmp$model$aicc,
AIC = hw_multipDmp$model$aic))
model_dlm = rbind(model_dlm,cbind(Model="hw_multiDmpExp",MASE=accuracy(hw_multiDmpExp)[6],
BIC = hw_multiDmpExp$model$bic,
AICC = hw_multiDmpExp$model$aicc,
AIC = hw_multiDmpExp$model$aic))
model_dlm = rbind(model_dlm,cbind(Model="hw_multipExp",MASE=accuracy(hw_multipExp)[6],
BIC = hw_multipExp$model$bic,
AICC = hw_multipExp$model$aicc,
AIC = hw_multipExp$model$aic))
model_dlm## Model MASE BIC AICC
## 1 hw_addi 0.247160046952369 5511.07634885595 5435.66154268864
## 2 hw_multip 0.223307745945573 6725.11394280828 6649.69913664097
## 3 hw_multipDmp 0.203561894680514 6447.56127723734 6367.76804289028
## 4 hw_multiDmpExp 0.203200910880705 6446.18950399393 6366.39626964687
## 5 hw_multipExp 0.23204037862927 6660.57583100187 6585.16102483456
## AIC
## 1 5434.7082716606
## 2 6648.74586561293
## 3 6366.70096020697
## 4 6365.32918696357
## 5 6584.20775380652hw_multipDmp = hw(solrad,seasonal="multiplicative" , damped = T)
checkresiduals(hw_multipDmp,test = F)
hw_multipExp = hw(solrad,seasonal="multiplicative" , damped = F, exponential = T)
checkresiduals(hw_multipExp,test = F)
ETS Methods
getEtsModelParam <- function(ts_obj){
err = c("N","A","M")
trnd = c("N","A","M")
sesn = c("N","A","M")
dmp = c(TRUE, FALSE)
mdl_mthd = expand.grid(err,trnd,sesn,dmp)
model_summ_df = data.frame(Model=character(), MASE=numeric(),
BIC=numeric(), AICC=numeric(), AIC = numeric())
for(i in 1:nrow(mdl_mthd)){
ets_mthd = paste(mdl_mthd[i,1],mdl_mthd[i,2],mdl_mthd[i,3],sep = "")
tryCatch({
if(ets_mthd == "NNN") {
ets_mdl = ets(ts_obj, damped = mdl_mthd[i,4], ic=c("aicc", "aic", "bic"))
model_summ_df = rbind(model_summ_df, cbind(Model=ets_mdl$method,
MASE=accuracy(ets_mdl)[6],
BIC = ets_mdl$bic,AICC = ets_mdl$aicc,
AIC = ets_mdl$aic))
} else if(mdl_mthd[i,2] != "N") {
ets_mdl = ets(ts_obj, model=ets_mthd,
damped=mdl_mthd[i,4], ic=c("aicc", "aic", "bic"))
model_summ_df = rbind(model_summ_df,
cbind(Model=ets_mdl$method,
MASE=accuracy(ets_mdl)[6],
BIC = ets_mdl$bic,AICC = ets_mdl$aicc,
AIC = ets_mdl$aic))
} else {
ets_mdl = ets(ts_obj, model=ets_mthd, damped=mdl_mthd[i,4],
ic=c("aicc", "aic", "bic"))
model_summ_df = rbind(model_summ_df, cbind(Model=ets_mdl$method,
MASE=accuracy(ets_mdl)[6],
BIC = ets_mdl$bic,AICC = ets_mdl$aicc,
AIC = ets_mdl$aic))
}
},error=function(e){
print(paste(mdl_mthd[i,],e,sep=":"))
})
}
return (model_summ_df)
}
getEtsModelParam(solrad)## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "TRUE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [1] "1:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Invalid error type\n"
## [1] "2:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [2] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [3] "3:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"
## [4] "FALSE:Error in ets(ts_obj, model = ets_mthd, damped = mdl_mthd[i, 4], ic = c(\"aicc\", : Forbidden model combination\n"## Model MASE BIC AICC
## 1 ETS(A,Ad,A) 0.246179687973231 5509.28207821646 5429.4888438694
## 2 ETS(A,Ad,N) 0.433469050524869 5959.4775686847 5932.6527667343
## 3 ETS(M,Ad,N) 0.658398691861859 6567.81938520252 6540.99458325212
## 4 ETS(M,Md,N) 0.66065437359335 6637.76980579802 6610.94500384762
## 5 ETS(A,Ad,A) 0.246179687973231 5509.28207821646 5429.4888438694
## 6 ETS(M,Ad,A) 0.37980945497033 6549.93960783857 6470.14637349151
## 7 ETS(M,Ad,M) 0.322257380231657 6034.36261969045 5954.56938534339
## 8 ETS(M,Md,M) 0.32011925393239 6076.41022774764 5996.61699340058
## 9 ETS(A,A,A) 0.247160046952369 5511.07634885595 5435.66154268864
## 10 ETS(A,N,N) 0.636820333080714 6309.84723126577 6296.40709712656
## 11 ETS(M,N,N) 0.636959906189964 6633.25257504203 6619.81244090282
## 12 ETS(A,A,N) 0.461151969887043 6026.25782306352 6003.88836700768
## 13 ETS(M,A,N) 0.651359673460806 6455.82498590279 6433.45552984695
## 14 ETS(M,M,N) 0.695815976506232 6631.62431864644 6609.2548625906
## 15 ETS(A,N,A) 0.2547040413186 5517.35758911071 5450.71933320031
## 16 ETS(M,N,A) 0.698942077110867 6614.23944723344 6547.60119132304
## 17 ETS(A,A,A) 0.247160046952369 5511.07634885595 5435.66154268864
## 18 ETS(M,A,A) 0.474856059737901 7679.12309973339 7603.70829356608
## 19 ETS(M,N,M) 0.313722575493707 6054.16199778708 5987.52374187668
## 20 ETS(M,A,M) 0.372166352645494 6182.32674943906 6106.91194327175
## 21 ETS(M,M,M) 0.529215103734309 6746.53608150663 6671.12127533931
## AIC
## 1 5428.42176118609
## 2 5932.52412967458
## 3 6540.8659461924
## 4 6610.8163667879
## 5 5428.42176118609
## 6 6469.0792908082
## 7 5953.50230266008
## 8 5995.54991071727
## 9 5434.7082716606
## 10 6296.37051176071
## 11 6619.77585553696
## 12 6003.79662388842
## 13 6433.36378672769
## 14 6609.16311947134
## 15 5449.9739915854
## 16 6546.85584970813
## 17 5434.7082716606
## 18 7602.75502253805
## 19 5986.77840026177
## 20 6105.95867224371
## 21 6670.16800431128etsauto <- ets(solrad,model="ZZZ")
etsauto## ETS(A,Ad,A)
##
## Call:
## ets(y = solrad, model = "ZZZ")
##
## 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.282ets = ets(solrad, damped = T, ic=c("aicc", "aic", "bic"), model = "AAA")
checkresiduals(ets, test=F)
dlm10 = dlm( y = as.vector(solrad) , x = as.vector(ppt), q = 10 )
polydlm10 = polyDlm( y = as.vector(solrad) , x = as.vector(ppt), q = 10, k= 2, show.beta = FALSE )
koyck= koyckDlm(y=as.vector(solrad),x=as.vector(ppt))
ardlm10 = ardlDlm( y = as.vector(solrad) , x = as.vector(ppt), p = 10, q= 10)
dynlm5 = dynlm(Y.t ~ L(Y.t ,1 ) + season(Y.t))
hw_multipExp = hw(solrad,seasonal="multiplicative" , damped = F, exponential = T)dynlmForecast = function(dynlm_model,y_series,x_series,frc_pd) {
q = frc_pd
Y.t = y_series
n = nrow(dynlm_model$model)+1
dynlm_model.frc = array(NA , (n + q))
dynlm_model.frc[1:n] = Y.t[1:length(Y.t)]
for (i in 1:q){
months = array(0,11)
months[(i-1)%%12] = 1
data.new = c(1,dynlm_model.frc[n-1+i],months)
dynlm_model.frc[n+i] = as.vector(dynlm_model$coefficients) %*% data.new
}
return(dynlm_model.frc)
}
dlmforecast = forecast(model = dlm10 , x = as.vector(forecastdata_ts) , h = 24)$forecasts
polyforecast = forecast(polydlm10 , x = as.vector(forecastdata_ts) , h = 24)$forecasts
koyckforecast = forecast(koyck, x = as.vector(forecastdata_ts), h = 24)$forecasts
ardlforecast = forecast(ardlm10, x = as.vector(forecastdata_ts) , h = 24)$forecasts
#dynlmforecast = forecast(dynlm5,solrad,ppt,24)
hwforecast = hw_multipExp$mean
etsforecast = as.vector(forecast(ets,h=24)$mean)
y_x = c(ppt,forecastdata_ts)
y_dlm = c(solrad,dlmforecast)
y_poly = c(solrad ,polyforecast)
y_koyck = c(solrad,koyckforecast)
y_ardl = c(solrad,ardlforecast)
#y_dynlm = frc_dynlm
y_hw = c(solrad, hwforecast)
y_ets = c(solrad,etsforecast)
Dataforecast = ts.intersect(
ts(y_x,start=c(1960,1),frequency = 12),
ts(y_dlm,start=c(1960,1),frequency = 12),
ts(y_poly,start=c(1960,1),frequency = 12),
ts(y_koyck,start=c(1960,1),frequency = 12),
ts(y_ardl,start=c(1960,1),frequency = 12),
#ts(y_dynlm,start=c(1960,1),frequency = 12),
ts(y_hw,start=c(1960,1),frequency = 12),
ts(y_ets,start=c(1960,1),frequency = 12)
)
colnames(Dataforecast) = c ("PPT", "DLM","Poly DLM", "Koyck", "ARDL", "Holt Winters" , "ETS")par(mar=c(4,4.5,5,1),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)
ts.plot(Dataforecast,xlim=c(2014,2016),ylim=c(0,40),
gpars= list(col=RColorBrewer::brewer.pal(7, "Set1")),
plot.type = c("single"),
main = "Forecasting of Solar Radtiation using the best fit of each model")
legend("topleft",text.width = 5, lty = 1, adj = c(0, 0.6),
bty = "n",col=RColorBrewer::brewer.pal(8, "Set1"),
cex = 0.65, c("PPT","DLM", "Polynomial", "Koyck", "ARDL", "HW" , "ETS"))
Task 2
task2 <- read.csv("data2.csv")
ppi <- ts(task2$price, start = c(2003,1), frequency = 4)
change <- ts(task2$change, start= c(2003,1), frequency = 4)
ppichange <- ts.intersect(ppi, change)
plot(ppichange, main= "Property Price v/s Change in population")
par(mar=c(4,4.5,5,1),cex.main=0.75, cex.lab=0.75, cex.axis=0.75)
ccf(as.vector(ppichange[,1]), as.vector(ppichange[,2]),
ylab='CCF', main = "Sample CCF between Residential Property Price Index and population change")
me.dif=ts.intersect(diff(diff(ppichange[,1],12)), diff(diff(ppichange[,2],12)))
prewhiten(as.vector(me.dif[,1]),as.vector(me.dif[,2]),ylab='CCF',
main="Sample CCF after prewhitening of the Residential Property Price Index and population change")