Task 1

Importing the file

data1 <- read_csv("data1.csv")

## Rows: 660 Columns: 2
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (2): solar, ppt
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

View(data1)

Converting the data set into Time series

data1 = ts(data1, start = c(1960, 1), end = c(2014,12), frequency = 12)
solar = data1[, 1]
ppt = data1[, 2]

Descriptive Anlaysis for Each time series for better understanding of the TS.

Stanarioty Check with ADF test

adf.test(solar)

## Warning in adf.test(solar): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  solar
## Dickey-Fuller = -4.7661, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

adf.test(ppt)

## Warning in adf.test(ppt): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ppt
## Dickey-Fuller = -6.7438, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

The results from ADF test concludes that both the time series, solar and ppt are non stationairty as the p-value is <0.05

Analysis

plot(solar,ylab='Solar Radition in W/m2',xlab='Year',
 type='o', main="Amount of Solar Radiation")

Seasonality is noticed in time series Solar Radiation

plot(ppt,ylab='Precipitation unit',xlab='Year',
 type='o', main="Precipitation Series")

Seasonality is noticed in time series Solar Radiation

Autocorrelation Check

acf(solar, main="ACF for Solar Radition in W/m2 ")

acf(ppt, main="ACF for Precipitation Unit")

For both the time series it its noted that significant number of lags are above the 95% of the confidence level, which proves non sationarity in both of the time series.

PACF Check

pacf(solar, main = "PACF for Solar Series",cex.main=0.05)

pacf(ppt, main = "PACF for PPT series",cex.main=0.05)

The first lag for both the series Solar and PPT is significant, which proves the series is non - stationrity.
After doing ADF, ACF and PACF testing we can conclude that our time series is non - stationrity.

Transformation for Solar TS

BC =  BoxCox.ar(solar)

BC$ci

## [1] 0.9 0.9

landings.tr = log (solar)

plot(landings.tr,ylab='Solar Radition in W/m2',xlab='Year',
 type='o', main="Amount of Solar Radition")

Seasonal Adjusting

par(mfrow=c(1,1))
landings.diff = diff(landings.tr, differences = 1, lag=12)
plot(landings.diff,ylab='Landings',
xlab='Year',main = "Time series plot of the first difference of monthly
landings.")
points(y=landings.diff,x=time(landings.diff), pch=as.vector(season(landings.diff)))

As seasonality was discovered in Solar time series, it was seasonally adjusted.

Transformation for Percipitation TS

BC =  BoxCox.ar(ppt)

BC$ci

## [1] 0.3 0.3

landings.tr2 = log (ppt)

plot(landings.tr,ylab='precipitation series',xlab='Year',
 type='o', main="Precipitation Unit")

Seasonal Adjusting

par(mfrow=c(1,1))
landings.diff2 = diff(landings.tr2, differences = 1, lag=12)
plot(landings.diff2,ylab='Landings Units',
xlab='Year',main = "Time series plot of the first difference of monthly
landings.")
points(y=landings.diff2,x=time(landings.diff2), pch=as.vector(season(landings.diff2)))

As seasonality was discovered in PPT series, it was seasonally adjusted.

Choice of Variable

For this forecasting, I have chosen PPT (X) as Independent Variable and Solar (Y) as Dependent Variable

selected_columns <- data1[, c(1, 2)]
data.scaled = scale(selected_columns)
plot (data.scaled, plot.type="s",col = c("blue", "red"),main = "Blue = Solar Radition, Red = Precipitation")

As there was a huge difference between means of solar Ts and PPT ts, I had to use “scaled data” function to show correlation

Testing of Models

Distributed Lag Models

model1 = dlm ( x = as.vector(ppt), y = as.vector (solar), q = 4)
summary(model1)

## 
## Call:
## lm(formula = model.formula, data = design)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.418  -5.743  -1.444   4.398  32.225 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.3727     0.7723  25.086  < 2e-16 ***
## x.t         -10.8482     1.7492  -6.202 9.93e-10 ***
## x.1           0.3550     2.5524   0.139    0.889    
## x.2          -0.2807     2.5946  -0.108    0.914    
## x.3          -0.4812     2.5406  -0.189    0.850    
## x.4           7.9026     1.7380   4.547 6.49e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.411 on 650 degrees of freedom
## Multiple R-squared:  0.2714, Adjusted R-squared:  0.2658 
## F-statistic: 48.42 on 5 and 650 DF,  p-value: < 2.2e-16
## 
## AIC and BIC values for the model:
##      AIC      BIC
## 1 4663.6 4695.003

Checking Residuals

checkresiduals(model1$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 586.43, df = 10, p-value < 2.2e-16

The ACF is way beyond the confidence value
The adjusted R Squared Value = 0.268
Only 1st and 4th lag is significiant

Polynomial Lag model

model1.2 = polyDlm(x = as.vector(ppt), y = as.vector (solar), q= 8, k = 2, show.beta = TRUE)

## Estimates and t-tests for beta coefficients:
##        Estimate Std. Error t value  P(>|t|)
## beta.0  -6.7500      0.574 -11.800 4.28e-29
## beta.1  -2.9000      0.344  -8.440 2.06e-16
## beta.2   0.0263      0.250   0.105 9.16e-01
## beta.3   2.0400      0.259   7.870 1.50e-14
## beta.4   3.1400      0.274  11.500 8.81e-28
## beta.5   3.3200      0.259  12.800 1.31e-33
## beta.6   2.5800      0.250  10.300 3.22e-23
## beta.7   0.9270      0.344   2.700 7.18e-03
## beta.8  -1.6400      0.574  -2.870 4.30e-03

summary (model1.2)

## 
## Call:
## "Y ~ (Intercept) + X.t"
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.818  -5.785  -1.357   4.376  31.750 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 17.43494    1.03121   16.91   <2e-16 ***
## z.t0        -6.75306    0.57384  -11.77   <2e-16 ***
## z.t1         4.30672    0.32482   13.26   <2e-16 ***
## z.t2        -0.45851    0.03978  -11.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.329 on 648 degrees of freedom
## Multiple R-squared:  0.2854, Adjusted R-squared:  0.2821 
## F-statistic: 86.25 on 3 and 648 DF,  p-value: < 2.2e-16

Checking Residuals

checkresiduals(model1.2$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 574.33, df = 10, p-value < 2.2e-16

For Model 1.2 the Z value for 1st, 2nd and 3rd lag are all significant
R.Sqaured value is 0.2821 which is good for model fitting
The ACF is beyond the confidence level, we can observe some seasonality

Koyck Model

model1.3 = koyckDlm(x = as.vector(ppt), y = as.vector (solar))
summary(model1.3)

## 
## Call:
## "Y ~ (Intercept) + Y.1 + X.t"
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.0926  -3.5961   0.3176   3.6103  14.8399 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.23925    0.76549  -2.925  0.00356 ** 
## Y.1          0.98546    0.02424  40.650  < 2e-16 ***
## X.t          5.34684    0.84383   6.336 4.37e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.814 on 656 degrees of freedom
## Multiple R-Squared: 0.7598,  Adjusted R-squared: 0.7591 
## Wald test:  1104 on 2 and 656 DF,  p-value: < 2.2e-16 
## 
## Diagnostic tests:
## NULL
## 
##                              alpha     beta       phi
## Geometric coefficients:  -154.0203 5.346844 0.9854613

Checking Residuals

checkresiduals(model1.3$model)

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 1413.2, df = 10, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 10

For model 1.3 we can interpret that the adjusted square value is better than other models i.e 0.7591
ACF values still over 95% confidence level

Autoregressive Lag Model

model1.4 = ardlDlm(x = as.vector(ppt), y = as.vector (solar), p= 1, q = 1)
summary(model1.4)

## 
## Time series regression with "ts" data:
## Start = 2, End = 660
## 
## Call:
## dynlm(formula = as.formula(model.text), data = data, start = 1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.6739  -2.8807  -0.3641   2.8687  20.1193 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.81174    0.53016   1.531    0.126    
## X.t         -6.99904    0.73480  -9.525   <2e-16 ***
## X.1          8.67630    0.71609  12.116   <2e-16 ***
## Y.1          0.91001    0.01851  49.161   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.027 on 655 degrees of freedom
## Multiple R-squared:  0.8321, Adjusted R-squared:  0.8314 
## F-statistic:  1082 on 3 and 655 DF,  p-value: < 2.2e-16

Checking Residuals

checkresiduals(model1.4$model)

## 
##  Breusch-Godfrey test for serial correlation of order up to 10
## 
## data:  Residuals
## LM test = 377.59, df = 10, p-value < 2.2e-16

In model 1.4 the Adjusted R - Squared value is most appropriate of all DLM models at 0.8314
No serious correlation
Residuals are not normally distributed, it has a slight right tail

Exponential Smoothing Models

model2 <- ses(solar, initial="simple", h=5)
summary(model2)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = solar, h = 5, initial = "simple") 
## 
##   Smoothing parameters:
##     alpha = 1 
## 
##   Initial states:
##     l = 5.0517 
## 
##   sigma:  4.5688
## Error measures:
##                        ME     RMSE      MAE       MPE     MAPE     MASE
## Training set 0.0001462894 4.568777 3.876091 -5.211851 27.29823 0.636771
##                   ACF1
## Training set 0.6677846
## 
## Forecasts:
##          Point Forecast      Lo 80    Hi 80      Lo 95    Hi 95
## Jan 2015        5.14828 -0.7068426 11.00340  -3.806357 14.10292
## Feb 2015        5.14828 -3.1321139 13.42867  -7.515490 17.81205
## Mar 2015        5.14828 -4.9930900 15.28965 -10.361607 20.65817
## Apr 2015        5.14828 -6.5619655 16.85853 -12.760995 23.05756
## May 2015        5.14828 -7.9441725 18.24073 -14.874898 25.17146

Checking Residuals

checkresiduals(model2)

## 
##  Ljung-Box test
## 
## data:  Residuals from Simple exponential smoothing
## Q* = 4227.3, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The data follows a normal distribution
The MASE value for model 2 is 0.636771

Holt Model

model2.1 <- holt(solar, initial="simple", h=5) 
summary(model2.1)

## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = solar, h = 5, initial = "simple") 
## 
##   Smoothing parameters:
##     alpha = 0.9165 
##     beta  = 1 
## 
##   Initial states:
##     l = 5.0517 
##     b = 1.3641 
## 
##   sigma:  3.6493
## Error measures:
##                        ME     RMSE      MAE      MPE     MAPE      MASE
## Training set -0.004941411 3.649286 2.806374 6.556498 21.13578 0.4610361
##                    ACF1
## Training set 0.06518525
## 
## Forecasts:
##          Point Forecast      Lo 80     Hi 80     Lo 95    Hi 95
## Jan 2015      3.3755382  -1.301210  8.052286  -3.77693 10.52801
## Feb 2015      1.7507183  -8.358924 11.860360 -13.71065 17.21208
## Mar 2015      0.1258983 -16.851847 17.103643 -25.83932 26.09112
## Apr 2015     -1.4989216 -26.473564 23.475721 -39.69434 36.69650
## May 2015     -3.1237415 -37.070990 30.823507 -55.04158 48.79409

Checking Residuals

checkresiduals(model2.1)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt's method
## Q* = 1096.3, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The MASE for model 2.1 is 0.4610361, which is relatively higher
The P value is less than <0.05 which means it’s significance.
Residuals is normally distributed
Acf is beyond the confidence level

Model with damped trend

model2.2 <- holt(solar, damped=TRUE, initial="simple", h=5) 
summary(model2.2)

## 
## Forecast method: Damped Holt's method
## 
## Model Information:
## Damped Holt's method 
## 
## Call:
##  holt(y = solar, h = 5, damped = TRUE, initial = "simple") 
## 
##   Smoothing parameters:
##     alpha = 0.9278 
##     beta  = 0.9278 
##     phi   = 0.8 
## 
##   Initial states:
##     l = 13.7194 
##     b = 2.5785 
## 
##   sigma:  3.4657
## 
##      AIC     AICc      BIC 
## 5932.524 5932.652 5959.477 
## 
## Error measures:
##                       ME     RMSE      MAE      MPE     MAPE     MASE
## Training set -0.00831309 3.452592 2.638613 3.791232 19.25658 0.433476
##                    ACF1
## Training set 0.07287524
## 
## Forecasts:
##          Point Forecast       Lo 80    Hi 80      Lo 95    Hi 95
## Jan 2015      3.7746889  -0.6668422  8.21622  -3.018047 10.56742
## Feb 2015      2.7358686  -5.9098357 11.38157 -10.486595 15.95833
## Mar 2015      1.9048121 -11.3560708 15.16570 -18.375958 22.18558
## Apr 2015      1.2399667 -16.7554651 19.23540 -26.281671 28.76160
## May 2015      0.7080903 -22.0017301 23.41791 -34.023583 35.43976

Checking Residuals

checkresiduals(model2.2)

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt's method
## Q* = 1134.4, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The MASE value for model2.2 is 0.433476
The P value is <0.05, which proves the non staionarity
Residuals are normally distributed.

Model with addictive Trend

model2.3 <- hw(solar,seasonal="additive", h=5*frequency(solar))
summary(model2.3)

## 
## Forecast method: Holt-Winters' additive method
## 
## Model Information:
## Holt-Winters' additive method 
## 
## Call:
##  hw(y = solar, h = 5 * frequency(solar), 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.5068095 10.74799
## Feb 2015       8.655602   4.3802222 12.930983   2.1169727 15.19423
## Mar 2015      14.122221   8.8814747 19.362968   6.1071910 22.13725
## Apr 2015      18.366625  12.3095941 24.423656   9.1031957 27.63005
## May 2015      23.522264  16.7439512 30.300578  13.1557290 33.88880
## Jun 2015      26.390991  18.9586817 33.823300  15.0242549 37.75773
## Jul 2015      27.219228  19.1837595 35.254697  14.9300392 39.50842
## Aug 2015      24.290519  15.6920179 32.889020  11.1402463 37.44079
## Sep 2015      19.495874  10.3670365 28.624711   5.5345219 33.45723
## Oct 2015      13.253765   3.6218809 22.885648  -1.4769303 27.98446
## Nov 2015       8.042144  -2.0695730 18.153861  -7.4223926 23.50668
## Dec 2015       5.840110  -4.7313921 16.411612 -10.3276072 22.00783
## Jan 2016       6.819616  -4.1942247 17.833457 -10.0246000 23.66383
## Feb 2016       9.347819  -2.0927719 20.788410  -8.1490551 26.84469
## Mar 2016      14.814438   2.9609172 26.667959  -3.3139578 32.94283
## Apr 2016      19.058842   6.8048110 31.312872   0.3179190 37.79976
## May 2016      24.214481  11.5711780 36.857784   4.8782175 43.55074
## Jun 2016      27.083208  14.0608586 40.105557   7.1672434 46.99917
## Jul 2016      27.911445  14.5194061 41.303485   7.4300887 48.39280
## Aug 2016      24.982736  11.2296053 38.735867   3.9491378 46.01633
## Sep 2016      20.188091   6.0818048 34.294377  -1.3856120 41.76179
## Oct 2016      13.945981  -0.5061083 28.398071  -8.1565824 36.04855
## Nov 2016       8.734361  -6.0566988 23.525420 -13.8866128 31.35533
## Dec 2016       6.532327  -8.5913309 21.655984 -16.5973115 29.66196
## Jan 2017       7.511833  -7.9385282 22.962194 -16.1174555 31.14112
## Feb 2017      10.040036  -5.7313777 25.811450 -14.0802599 34.16033
## Mar 2017      15.506655  -0.5805619 31.593872  -9.0966201 40.10993
## Apr 2017      19.751059   3.3529822 36.149135  -5.3276351 44.82975
## May 2017      24.906698   8.2024281 41.610968  -0.6402783 50.45367
## Jun 2017      27.775425  10.7693726 44.781476   1.7669125 53.78394
## Jul 2017      28.603662  11.3000071 45.907317   2.1400055 55.06732
## Aug 2017      25.674953   8.0776596 43.272246  -1.2377847 52.58769
## Sep 2017      20.880308   2.9931440 38.767471  -6.4757484 48.23636
## Oct 2017      14.638198  -3.5352503 32.811647 -13.1556929 42.43209
## Nov 2017       9.426578  -9.0297392 27.882894 -18.7999231 37.65308
## Dec 2017       7.224543 -11.5113813 25.960468 -21.4295807 35.87867
## Jan 2018       8.204050 -10.8084216 27.216522 -20.8730161 37.28112
## Feb 2018      10.732253  -8.5537326 30.018238 -18.7631166 40.22762
## Mar 2018      16.198872  -3.3577754 35.755519 -13.7104391 46.10818
## Apr 2018      20.443275   0.6186997 40.267851  -9.8757968 50.76235
## May 2018      25.598915   5.5090332 45.688797  -5.1259078 56.32374
## Jun 2018      28.467641   8.1149718 48.820311  -2.6590806 59.59436
## Jul 2018      29.295879   8.6828417 49.908916  -2.2290411 60.82080
## Aug 2018      26.367170   5.4960925 47.238247  -5.5523883 58.28673
## Sep 2018      21.572524   0.4456486 42.699400 -10.7382439 53.88329
## Oct 2018      15.330415  -6.0501007 36.710931 -17.3682620 48.02909
## Nov 2018      10.118794 -11.5132798 31.750868 -22.9646081 43.20220
## Dec 2018       7.916760 -13.9648644 29.798385 -25.5482967 41.38182
## Jan 2019       8.896267 -13.2330166 31.025550 -24.9475516 42.74009
## Feb 2019      11.424470 -10.9505524 33.799492 -22.7951736 45.64411
## Mar 2019      16.891089  -5.7278620 39.510039 -17.7016112 51.48379
## Apr 2019      21.135492  -1.7256366 43.996621 -13.8275871 56.09857
## May 2019      26.291132   3.1895190 49.392744  -9.0397360 61.62200
## Jun 2019      29.159858   5.8194018 52.500314  -6.5362894 64.85601
## Jul 2019      29.988096   6.4103847 53.565807  -6.0709017 66.04709
## Aug 2019      27.059386   3.2459605 50.872812  -9.3601057 63.47888
## Sep 2019      22.264741  -1.7829062 46.312389 -14.5129618 59.04244
## Oct 2019      16.022632  -8.2577885 40.303052 -21.1110666 53.15633
## Nov 2019      10.811011 -13.7007762 35.322798 -26.6765326 48.29855
## Dec 2019       8.608977 -16.1328122 33.350766 -29.2303243 46.44828

Checking Residuals

checkresiduals(model2.3)

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' additive method
## Q* = 205.55, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The Mase Value for model2.3 is 0.24716
The Residuals shows a right tail which means it is not a normal distribution

State Space Models

model3  = ets(solar, model="ANN")
summary(model3)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = solar, 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

Checking Residuals

checkresiduals(model3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,N,N)
## Q* = 4227.6, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

MASE for model3 is 0.6368203
Residuals is normally distributed
ACF is beyond confidence level

model3.1 = ets(solar, model="MNN")
summary(model3.1)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = solar, 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.6678785

Checking Residuals

checkresiduals(model3.1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,N,N)
## Q* = 1334.8, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

MASE for model 3.1 is 0.6369599
We can see a right tail in residuals, hence it is not normally distributed

model3.2 = ets(solar, model="AAN")
summary(model3.2)

## ETS(A,Ad,N) 
## 
## Call:
##  ets(y = solar, 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

Checking Residuals

checkresiduals(model3.2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,N)
## Q* = 1134, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

MASE for model 3.2 is 0.4334691, which is an improvement
Residuals is normally distributed
P - value is less than 0.05, which shows correlation
ACF values are beyond the confidence level

model3.3 = ets(solar, model="MAN", damped = TRUE)
summary(model3.3)

## ETS(M,Ad,N) 
## 
## Call:
##  ets(y = solar, model = "MAN", damped = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9815 
##     beta  = 1e-04 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 9.9304 
##     b = 5.8356 
## 
##   sigma:  0.35
## 
##      AIC     AICc      BIC 
## 6540.866 6540.995 6567.819 
## 
## Training set error measures:
##                      ME    RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set -0.4464292 4.75387 4.007741 -9.78188 29.60799 0.6583987 0.6870915

Checking Residuals

checkresiduals(model3.3)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Ad,N)
## Q* = 1745.9, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

MASE value for model 3.3 is 0.6583987
Residuals has a right tail, which means it is not normally distributed

model3.4 = ets(solar, model="AAA")
summary(model3.4)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = solar, 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.1700724

Checking Residuals

checkresiduals(model3.4)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 210.76, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The MASE for model 3.4 0.2461797, which is the most appropriate model by now
The P value is <0.5, which proves non - stationarity
The ACF value is mostly adjusted
The residuals is normally distributed.

model3.5 = ets(solar,model = "MMN")
summary(model3.5)

## ETS(M,M,N) 
## 
## Call:
##  ets(y = solar, 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

Checking Residuals

checkresiduals(model3.5)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,M,N)
## Q* = 1371.3, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

MASE for model 3.5 is 0.695816
ACF is beyond confidence level
Residuals has a right tail, which means it is not normally distributted

model3.6 = ets(solar,model = "MMM")
summary(model3.6)

## ETS(M,Md,M) 
## 
## Call:
##  ets(y = solar, model = "MMM") 
## 
##   Smoothing parameters:
##     alpha = 0.7892 
##     beta  = 0.016 
##     gamma = 0.0696 
##     phi   = 0.9575 
## 
##   Initial states:
##     l = 9.2782 
##     b = 1.1745 
##     s = 0.7133 0.3063 0.6852 1.1282 1.2133 1.5082
##            1.7211 1.2962 1.1358 0.9309 0.656 0.7056
## 
##   sigma:  0.2341
## 
##      AIC     AICc      BIC 
## 5995.550 5996.617 6076.410 
## 
## Training set error measures:
##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.02448728 3.076938 1.948599 -5.645612 16.87559 0.3201193
##                   ACF1
## Training set 0.1452617

Checking Residuals

checkresiduals(model3.6)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,Md,M)
## Q* = 79.153, df = 24, p-value = 8.298e-08
## 
## Model df: 0.   Total lags used: 24

MASE value for model3.6 is 0.3201193
The residuals are not normally distributed

Calculating MASE Values for Time series regression (DLM Models)

MASE(model1, model1.2, model1.3, model1.4)

##            n      MASE
## model1   656 1.6463575
## model1.2 652 1.6190338
## model1.3 659 1.0324829
## model1.4 659 0.8392434

Implementation of model for Forecasting two years of Solar Radiation

After considering MASE values for every possible model. It comes down to conclusion to forecast on model 3.4 has it has the lowest MASE value of all
model3.4 is a Local Additive Seasonal Model (AAA) which will be used for forecasting

solarfor =forecast::forecast(model3.4, h = 24)
summary(solarfor)

## 
## Forecast method: ETS(A,Ad,A)
## 
## Model Information:
## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = solar, 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 
## 
## 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

plot(solarfor)

A visual presentation of Solar forecasting

Forecasting of Solar Radition till 2018

#Importing the CSV files for PPT values

datax <- read_csv("data.x.csv")

## Rows: 24 Columns: 1
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (1): x
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Taking out the forecasted values of Solar till 2016

solar.val = c(solarfor$mean)

#Renaming the column in datax as PPT for better understanding

colnames(datax) = "PPT"

Making a dataset of Forecasted valuesof Solar and PPT values till 2016

ogdata = data.frame(solar = solar.val, ppt = datax)

Converting the data set into Time Series

data.x = ts(ogdata, start = c(2015, 1), end = c(2016,12), frequency = 12)
solar1 = data1[, 1]
ppt1 = data1[, 2]

Stationarity Check

adf.test(solar1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  solar1
## Dickey-Fuller = -4.7661, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

adf.test(ppt1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ppt1
## Dickey-Fuller = -6.7438, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary

As per ADF test the p - value is smaller than 0.05 which means our new time series is non stationary

Model Checking - As dicvoered earlier Local Additive Seasonal Model is best suited model for our forecasting

modela = ets(solar1, model="AAA")
summary(modela)

## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = solar1, 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.1700724

Check Residuals

checkresiduals(modela)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,A)
## Q* = 210.76, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

The MASE value for the selected model is 0.2461797 which was the lowest as per my comparison
The ACF shows correlation
The Residuals is normally distributed

Forecasting of Solar Radition till 2018 on our selected model

finalfor <- forecast::forecast(modela, h = 24 * 2)
summary(finalfor)

## 
## Forecast method: ETS(A,Ad,A)
## 
## Model Information:
## ETS(A,Ad,A) 
## 
## Call:
##  ets(y = solar1, 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 
## 
## 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.9405301  8.949866   1.3499548 10.54044
## Feb 2015       8.479778   4.2305479 12.729009   1.9811413 14.97842
## Mar 2015      13.784309   8.5799016 18.988716   5.8248548 21.74376
## Apr 2015      17.598568  11.5887764 23.608360   8.4073847 26.78975
## May 2015      22.632258  15.9128030 29.351713  12.3557383 32.90878
## Jun 2015      25.599150  18.2380236 32.960277  14.3412786 36.85702
## Jul 2015      26.761775  18.8104968 34.713053  14.6013444 38.92221
## Aug 2015      23.722251  15.2216103 32.222892  10.7216429 36.72286
## Sep 2015      18.222624   9.2059552 27.239293   4.4328191 32.01243
## Oct 2015      12.293259   2.7884705 21.798047  -2.2430606 26.82958
## Nov 2015       7.504219  -2.4648770 17.473315  -7.7421977 22.75064
## Dec 2015       5.150488  -5.2622860 15.563263 -10.7744758 21.07545
## Jan 2016       5.947275  -4.8911752 16.785726 -10.6287044 22.52326
## Feb 2016       8.481728  -2.7662508 19.729707  -8.7205712 25.68403
## Mar 2016      13.786139   2.1429879 25.429291  -4.0205242 31.59280
## Apr 2016      17.600287   5.5749058 29.625668  -0.7909464 35.99152
## May 2016      22.633871  10.2380087 35.029734   3.6760353 41.59171
## Jun 2016      25.600665  12.8450464 38.356283   6.0926298 45.10870
## Jul 2016      26.763197  13.6576670 39.868726   6.7200186 46.80637
## Aug 2016      23.723586  10.2772227 37.169950   3.1591478 44.28802
## Sep 2016      18.223877   4.4450855 32.002669  -2.8489663 39.29672
## Oct 2016      12.294435  -1.8089722 26.397843  -9.2748652 33.86374
## Nov 2016       7.505324  -6.9154138 21.926061 -14.5492911 29.55994
## Dec 2016       5.151525  -9.5797258 19.882776 -17.3779790 27.68103
## Jan 2017       5.948249  -9.0871906 20.983688 -17.0464714 28.94297
## Feb 2017       8.482642  -6.8508988 23.816183 -14.9679850 31.93327
## Mar 2017      13.786997  -1.8389729 29.412968 -10.1108619 37.68486
## Apr 2017      17.601092   1.6880528 33.514132  -6.7358014 41.93799
## May 2017      22.634628   6.4395949 38.829660  -2.1335375 47.40279
## Jun 2017      25.601375   9.1291648 42.073585   0.4093036 50.79345
## Jul 2017      26.763863  10.0190534 43.508673   1.1548866 52.37284
## Aug 2017      23.724212   6.7111602 40.737263  -2.2950052 49.74343
## Sep 2017      18.224465   0.9473269 35.501602  -8.1986374 44.64757
## Oct 2017      12.294987  -5.2422687 29.832242 -14.5259309 39.11590
## Nov 2017       7.505841 -10.2877372 25.299420 -19.7070886 34.71877
## Dec 2017       5.152011 -12.8942567 23.198279 -22.4473739 32.75140
## Jan 2018       5.948705 -12.3468264 24.244236 -22.0318957 33.92931
## Feb 2018       8.483070 -10.0583229 27.024464 -19.8735437 36.83968
## Mar 2018      13.787399  -4.9966434 32.571442 -14.9403150 42.51511
## Apr 2018      17.601470  -1.4221329 36.625072 -11.4926198 46.69556
## May 2018      22.634982   3.3747944 41.895169  -6.8209329 52.09090
## Jun 2018      25.601707   6.1078017 45.095613  -4.2116486 55.41506
## Jul 2018      26.764175   7.0393167 46.489034  -3.4023928 56.93074
## Aug 2018      23.724505   3.7713624 43.677647  -6.7911931 54.24020
## Sep 2018      18.224740  -1.9541074 38.403587 -12.6361439 49.08562
## Oct 2018      12.295245  -8.1068132 32.697304 -18.9070105 43.49750
## Nov 2018       7.506084 -13.1167729 28.128941 -24.0338538 39.04602
## Dec 2018       5.152239 -15.6890799 25.993558 -26.7218076 37.02629

Forecasting for Jan 2016 to Dec 2018 were made

Visual Presenation for Forecasting

plot(finalfor)

The blue graphical area shows the forecasting of Solar Radiations.

Task 2

#Importing the Dataset

data2 <- read_csv("data2.csv")

## Rows: 54 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Quarter
## dbl (2): price, change
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

str(data2)

## spc_tbl_ [54 × 3] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Quarter: chr [1:54] "Sep-2003" "Dec-2003" "Mar-2004" "Jun-2004" ...
##  $ price  : num [1:54] 60.7 62.1 60.8 60.9 60.9 62.4 62.5 63.2 63.1 64 ...
##  $ change : num [1:54] 14017 12350 17894 9079 16210 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Quarter = col_character(),
##   ..   price = col_double(),
##   ..   change = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

Converting the data set into Time series

price =ts(data2$price, start = c(2003, 9), frequency = 4)
change = ts(data2$change, start = c(2003, 9), frequency = 4)

Establishing Coorelation Between Price and Change

change.joint=ts.intersect(price,change)
plot(change.joint,yax.flip=T)

Sample CFF between PPI and population change Time Series

ccf(as.vector(change.joint[,1]), as.vector(change.joint[,2]),ylab='CCF',
main = "Sample CCF between PPI and Population Change in Melbourne")

We observe substantial values that exceed the 95% confidence threshold, indicating the presence of a spurious correlation.

Prewhitening

pre.white=ts.intersect(diff(diff(price,4)),diff(diff(change,4)))
prewhiten(as.vector(pre.white [,1]),as.vector(pre.white[,2]),ylab='CCF',
main="Sample CFF after prewhitening ")

Before doing pre whitening, differencing is applied at 4 lags to remove any stationarity, if any and prewhitening is done on the difference model.
After pre whitening the spurious correlation has been removed.

Assignment 2

Geetesh Sanjeev Matreja (s3931386)

2023-09-20

Task 1

Importing the file

Converting the data set into Time series

Descriptive Anlaysis for Each time series for better understanding of the TS.

Stanarioty Check with ADF test

Analysis

Autocorrelation Check

PACF Check

Transformation for Solar TS

Seasonal Adjusting

Transformation for Percipitation TS

Seasonal Adjusting

Choice of Variable

Testing of Models

Distributed Lag Models

Checking Residuals

Polynomial Lag model

Checking Residuals

Koyck Model

Checking Residuals

Autoregressive Lag Model

Checking Residuals

Exponential Smoothing Models

Checking Residuals

Holt Model

Checking Residuals

Model with damped trend

Checking Residuals

Model with addictive Trend

Checking Residuals

State Space Models

Checking Residuals

Checking Residuals

Checking Residuals

Checking Residuals

Checking Residuals

Checking Residuals

Checking Residuals

Calculating MASE Values for Time series regression (DLM Models)

Implementation of model for Forecasting two years of Solar Radiation

Forecasting of Solar Radition till 2018

Taking out the forecasted values of Solar till 2016

Making a dataset of Forecasted valuesof Solar and PPT values till 2016

Converting the data set into Time Series

Stationarity Check

Model Checking - As dicvoered earlier Local Additive Seasonal Model is best suited model for our forecasting

Check Residuals

Forecasting of Solar Radition till 2018 on our selected model

Visual Presenation for Forecasting

Task 2

Converting the data set into Time series

Establishing Coorelation Between Price and Change

Sample CFF between PPI and population change Time Series

Prewhitening