##Libraries used
library(fpp2)
## Loading required package: ggplot2
## Loading required package: forecast
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Loading required package: fma
## Loading required package: expsmooth
library(ggplot2)
library(forecast)
library(fGarch)
## Loading required package: timeDate
## Loading required package: timeSeries
## Loading required package: fBasics
library(rugarch)
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
library(tseries)
#Data set
str(qcement)
##  Time-Series [1:233] from 1956 to 2014: 0.465 0.532 0.561 0.57 0.529 0.604 0.603 0.582 0.554 0.62 ...
head(qcement)
##       Qtr1  Qtr2  Qtr3  Qtr4
## 1956 0.465 0.532 0.561 0.570
## 1957 0.529 0.604

SEASONAL ARIMA MODEL

##Setting lag of 4 to represent quarters of production
autoplot(qcement)+ylab("Total quarterly production (in tonnes)")+xlab("1956:Q1 to 2014:Q1.")

qcement %>% diff(lag=4) %>% ggtsdisplay()

qcement %>% diff(lag=4) %>% diff() %>% ggtsdisplay()

qcement %>%
  Arima(order = c(0,1,1), seasonal = c(0,1,1)) %>% residuals() %>% ggtsdisplay()

##
A1 <- Arima(qcement, order=c(0,1,3), seasonal=c(0,1,1))
checkresiduals(A1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,3)(0,1,1)[4]
## Q* = 12.54, df = 4, p-value = 0.01375
## 
## Model df: 4.   Total lags used: 8
A1 %>% forecast(h=4) %>% autoplot()

ETS MODEL

qcement.ets <- ets(qcement)
qcement
##       Qtr1  Qtr2  Qtr3  Qtr4
## 1956 0.465 0.532 0.561 0.570
## 1957 0.529 0.604 0.603 0.582
## 1958 0.554 0.620 0.646 0.637
## 1959 0.573 0.673 0.690 0.681
## 1960 0.621 0.698 0.753 0.728
## 1961 0.688 0.737 0.742 0.692
## 1962 0.637 0.757 0.783 0.757
## 1963 0.674 0.774 0.835 0.838
## 1964 0.797 0.904 0.949 0.975
## 1965 0.902 0.974 0.969 0.967
## 1966 0.849 0.961 0.966 0.922
## 1967 0.836 0.998 1.025 0.971
## 1968 0.892 0.973 1.047 1.017
## 1969 0.948 1.032 1.190 1.136
## 1970 1.049 1.134 1.229 1.188
## 1971 1.058 1.209 1.199 1.253
## 1972 1.070 1.282 1.303 1.281
## 1973 1.148 1.305 1.342 1.452
## 1974 1.184 1.352 1.316 1.353
## 1975 1.121 1.297 1.318 1.281
## 1976 1.109 1.299 1.341 1.290
## 1977 1.101 1.284 1.321 1.317
## 1978 1.122 1.261 1.312 1.298
## 1979 1.205 1.302 1.377 1.359
## 1980 1.232 1.386 1.440 1.439
## 1981 1.282 1.573 1.533 1.641
## 1982 1.337 1.575 1.475 1.357
## 1983 1.086 1.158 1.279 1.313
## 1984 1.166 1.373 1.456 1.495
## 1985 1.251 1.456 1.631 1.555
## 1986 1.375 1.546 1.568 1.561
## 1987 1.332 1.458 1.501 1.615
## 1988 1.418 1.625 1.770 1.791
## 1989 1.621 1.719 1.972 1.894
## 1990 1.565 1.645 1.658 1.668
## 1991 1.343 1.441 1.444 1.497
## 1992 1.289 1.501 1.539 1.568
## 1993 1.450 1.668 1.648 1.863
## 1994 1.468 1.755 1.962 1.833
## 1995 1.626 1.703 1.733 1.545
## 1996 1.526 1.593 1.706 1.699
## 1997 1.511 1.785 1.826 1.830
## 1998 1.719 1.861 1.956 2.067
## 1999 1.737 1.944 2.005 2.027
## 2000 1.835 2.070 1.898 1.652
## 2001 1.554 1.717 1.679 1.836
## 2002 1.729 1.992 2.030 1.978
## 2003 1.831 1.892 2.227 2.090
## 2004 1.963 2.180 2.307 2.157
## 2005 1.980 2.481 2.340 2.265
## 2006 2.027 2.278 2.427 2.451
## 2007 2.140 2.362 2.536 2.562
## 2008 2.183 2.558 2.612 2.373
## 2009 1.963 2.160 2.325 2.273
## 2010 1.904 2.401 2.494 2.296
## 2011 2.055 2.273 2.499 2.390
## 2012 2.067 2.223 2.451 2.503
## 2013 2.049 2.528 2.637 2.565
## 2014 2.229
checkresiduals(qcement.ets)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,A,M)
## Q* = 25.376, df = 3, p-value = 1.289e-05
## 
## Model df: 8.   Total lags used: 11
auto.arima(qcement)
## Series: qcement 
## ARIMA(2,0,2)(0,1,1)[4] with drift 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     sma1   drift
##       1.6175  -0.7194  -1.0358  0.3658  -0.7690  0.0082
## s.e.  0.1504   0.1236   0.1562  0.0755   0.0439  0.0011
## 
## sigma^2 estimated as 0.006588:  log likelihood=251.1
## AIC=-488.2   AICc=-487.69   BIC=-464.16
#Forecasting
E1 <- qcement %>% ets() %>% forecast(h=4) %>% autoplot()
E1

GARCH MODEL

qcement.garch <- garch(qcement)
## 
##  ***** ESTIMATION WITH ANALYTICAL GRADIENT ***** 
## 
## 
##      I     INITIAL X(I)        D(I)
## 
##      1     2.620279e-01     1.000e+00
##      2     5.000000e-02     1.000e+00
##      3     5.000000e-02     1.000e+00
## 
##     IT   NF      F         RELDF    PRELDF    RELDX   STPPAR   D*STEP   NPRELDF
##      0    1  5.457e+02
##      1    2  2.069e+02  6.21e-01  7.67e+00  9.0e-01  4.2e+03  1.0e+00  1.61e+04
##      2    4  2.040e+02  1.42e-02  9.82e-03  2.3e-02  2.0e+00  5.0e-02  5.78e+00
##      3    6  1.993e+02  2.29e-02  2.28e-02  4.1e-02  2.0e+00  1.0e-01  3.57e-01
##      4    8  1.986e+02  3.58e-03  3.59e-03  8.0e-03  1.2e+01  2.0e-02  2.59e-01
##      5   10  1.974e+02  6.08e-03  6.10e-03  1.6e-02  2.2e+00  4.0e-02  2.09e-01
##      6   12  1.972e+02  1.05e-03  1.05e-03  3.1e-03  2.0e+01  8.0e-03  1.83e-01
##      7   14  1.971e+02  2.02e-04  2.02e-04  6.1e-04  8.8e+01  1.6e-03  1.71e-01
##      8   17  1.968e+02  1.53e-03  1.53e-03  4.9e-03  3.7e+00  1.3e-02  1.68e-01
##      9   20  1.968e+02  2.92e-05  2.92e-05  9.6e-05  4.9e+02  2.6e-04  1.64e-01
##     10   22  1.968e+02  5.82e-05  5.82e-05  1.9e-04  6.2e+01  5.1e-04  1.61e-01
##     11   24  1.968e+02  1.16e-05  1.16e-05  3.8e-05  1.2e+03  1.0e-04  1.61e-01
##     12   26  1.968e+02  2.32e-06  2.32e-06  7.7e-06  6.1e+03  2.0e-05  1.60e-01
##     13   28  1.968e+02  4.64e-06  4.64e-06  1.5e-05  7.6e+02  4.1e-05  1.60e-01
##     14   30  1.968e+02  9.28e-07  9.28e-07  3.1e-06  1.5e+04  8.2e-06  1.60e-01
##     15   33  1.968e+02  1.86e-08  1.86e-08  6.1e-08  7.6e+05  1.6e-07  1.60e-01
##     16   35  1.968e+02  3.71e-08  3.71e-08  1.2e-07  9.5e+04  3.3e-07  1.60e-01
##     17   37  1.968e+02  7.43e-08  7.43e-08  2.5e-07  4.7e+04  6.6e-07  1.60e-01
##     18   39  1.968e+02  1.49e-08  1.49e-08  4.9e-08  9.5e+05  1.3e-07  1.60e-01
##     19   41  1.968e+02  2.97e-09  2.97e-09  9.8e-09  4.7e+06  2.6e-08  1.60e-01
##     20   43  1.968e+02  5.94e-09  5.94e-09  2.0e-08  5.9e+05  5.2e-08  1.60e-01
##     21   45  1.968e+02  1.19e-09  1.19e-09  3.9e-09  1.2e+07  1.0e-08  1.60e-01
##     22   47  1.968e+02  2.38e-10  2.38e-10  7.9e-10  5.9e+07  2.1e-09  1.60e-01
##     23   50  1.968e+02  1.90e-09  1.90e-09  6.3e-09  1.9e+06  1.7e-08  1.60e-01
##     24   53  1.968e+02  3.80e-11  3.80e-11  1.3e-10  1.5e+00  3.4e-10 -4.20e-02
##     25   55  1.968e+02  7.61e-11  7.61e-11  2.5e-10  1.5e+00  6.7e-10 -4.20e-02
##     26   58  1.968e+02  1.52e-12  1.52e-12  5.0e-12  1.5e+00  1.3e-11 -4.20e-02
##     27   61  1.968e+02  1.22e-11  1.22e-11  4.0e-11  1.5e+00  1.1e-10 -4.20e-02
##     28   65  1.968e+02  2.40e-14  2.43e-14  8.1e-14  1.5e+00  2.1e-13 -4.20e-02
##     29   67  1.968e+02  4.91e-14  4.87e-14  1.6e-13  1.5e+00  4.3e-13 -4.20e-02
##     30   69  1.968e+02  9.82e-15  9.73e-15  3.2e-14  1.5e+00  8.6e-14 -4.20e-02
##     31   71  1.968e+02  1.89e-14  1.95e-14  6.4e-14  1.5e+00  1.7e-13 -4.18e-02
##     32   73  1.968e+02  4.48e-15  3.89e-15  1.3e-14  1.5e+00  3.4e-14 -4.20e-02
##     33   76  1.968e+02  2.57e-14  2.65e-14  8.8e-14  1.5e+00  2.3e-13 -4.18e-02
##     34   78  1.968e+02 -5.08e+07  5.30e-15  1.8e-14  1.5e+00  4.7e-14 -4.20e-02
## 
##  ***** FALSE CONVERGENCE *****
## 
##  FUNCTION     1.968002e+02   RELDX        1.753e-14
##  FUNC. EVALS      78         GRAD. EVALS      34
##  PRELDF       5.299e-15      NPRELDF     -4.196e-02
## 
##      I      FINAL X(I)        D(I)          G(I)
## 
##      1    4.751213e-01     1.000e+00     1.335e+01
##      2    8.913371e-01     1.000e+00     9.858e+00
##      3    3.017120e-14     1.000e+00     1.490e+01
qcement.garch
## 
## Call:
## garch(x = qcement)
## 
## Coefficient(s):
##        a0         a1         b1  
## 4.751e-01  8.913e-01  3.017e-14
plot(qcement.garch)

summary(qcement.garch)
## 
## Call:
## garch(x = qcement)
## 
## Model:
## GARCH(1,1)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## 0.6049 0.8362 0.9368 1.0167 1.2471 
## 
## Coefficient(s):
##     Estimate  Std. Error  t value Pr(>|t|)
## a0 4.751e-01   1.085e+00    0.438    0.662
## a1 8.913e-01   1.412e+00    0.631    0.528
## b1 3.017e-14   1.240e+00    0.000    1.000
## 
## Diagnostic Tests:
##  Jarque Bera Test
## 
## data:  Residuals
## X-squared = 2.4635, df = 2, p-value = 0.2918
## 
## 
##  Box-Ljung test
## 
## data:  Squared.Residuals
## X-squared = 32.034, df = 1, p-value = 1.515e-08
qcement.fit <- garchFit(formula=~garch(1,1),qcement)
## 
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          norm
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          233
##  Recursion Init:            mci
##  Series Scale:              0.5395759
## 
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                       U         V   params includes
##     mu     -27.47285739  27.47286 2.747286     TRUE
##     omega    0.00000100 100.00000 0.100000     TRUE
##     alpha1   0.00000001   1.00000 0.100000     TRUE
##     gamma1  -0.99999999   1.00000 0.100000    FALSE
##     beta1    0.00000001   1.00000 0.800000     TRUE
##     delta    0.00000000   2.00000 2.000000    FALSE
##     skew     0.10000000  10.00000 1.000000    FALSE
##     shape    1.00000000  10.00000 4.000000    FALSE
##  Index List of Parameters to be Optimized:
##     mu  omega alpha1  beta1 
##      1      2      3      5 
##  Persistence:                  0.9 
## 
## 
## --- START OF TRACE ---
## Selected Algorithm: nlminb 
## 
## R coded nlminb Solver: 
## 
##   0:     290.24840:  2.74729 0.100000 0.100000 0.800000
##   1:     270.96138:  2.71736 0.0175098 0.111846 0.760830
##   2:     259.43079:  2.64869 0.0324435 0.199188 0.771846
##   3:     251.61092:  2.64108 0.00929393 0.197331 0.765766
##   4:     250.89873:  2.58379 0.0146697 0.207599 0.769620
##   5:     249.17855:  2.53506 0.00160820 0.217073 0.772623
##   6:     247.66068:  2.47887 0.0106005 0.218923 0.763602
##   7:     246.06071:  2.42763 0.00713435 0.229172 0.752693
##   8:     245.80748:  2.42778 0.00575023 0.229242 0.752599
##   9:     245.75732:  2.42745 0.00468039 0.230121 0.752610
##  10:     245.71799:  2.42686 0.00506603 0.231438 0.752656
##  11:     245.61875:  2.42513 0.00436533 0.237241 0.752252
##  12:     245.56622:  2.42785 0.00538995 0.241103 0.748035
##  13:     245.49342:  2.43092 0.00504065 0.244536 0.743391
##  14:     244.83871:  2.42809 0.00855788 0.304654 0.673637
##  15:     244.70387:  2.41026 0.00691775 0.310824 0.672640
##  16:     244.65516:  2.41698 0.00724173 0.319646 0.672641
##  17:     244.64739:  2.41692 0.00687046 0.319656 0.672581
##  18:     244.63995:  2.41367 0.00670697 0.320196 0.671839
##  19:     244.62950:  2.41345 0.00703059 0.322614 0.670069
##  20:     244.55031:  2.41850 0.00753002 0.346466 0.650363
##  21:     244.48548:  2.41685 0.00874040 0.364119 0.624918
##  22:     244.46823:  2.40291 0.00867575 0.388833 0.606909
##  23:     244.45743:  2.40923 0.00909669 0.387052 0.605304
##  24:     244.45738:  2.40915 0.00913806 0.387202 0.604982
##  25:     244.45737:  2.40914 0.00915181 0.387329 0.604851
##  26:     244.45737:  2.40914 0.00915280 0.387349 0.604843
##  27:     244.45737:  2.40914 0.00915276 0.387351 0.604844
## 
## Final Estimate of the Negative LLH:
##  LLH:  100.7029    norm LLH:  0.4322014 
##          mu       omega      alpha1       beta1 
## 1.299915995 0.002664753 0.387350565 0.604843605 
## 
## R-optimhess Difference Approximated Hessian Matrix:
##                 mu       omega      alpha1       beta1
## mu     -1910.76475   18032.631    21.44741    269.8891
## omega  18032.63088 -896627.013 -8251.84891 -18596.6710
## alpha1    21.44741   -8251.849  -613.31032   -675.3555
## beta1    269.88912  -18596.671  -675.35545   -928.4478
## attr(,"time")
## Time difference of 0.004199028 secs
## 
## --- END OF TRACE ---
## 
## 
## Time to Estimate Parameters:
##  Time difference of 0.0236001 secs
qcement.fit
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = qcement) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7f8470acd5c8>
##  [data = qcement]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##        mu      omega     alpha1      beta1  
## 1.2999160  0.0026648  0.3873506  0.6048436  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## mu      1.299916    0.025863   50.262  < 2e-16 ***
## omega   0.002665    0.001890    1.410  0.15861    
## alpha1  0.387351    0.121749    3.182  0.00146 ** 
## beta1   0.604844    0.120393    5.024 5.06e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -100.7029    normalized:  -0.4322014 
## 
## Description:
##  Wed Jul 29 19:34:56 2020 by user:
### plot(qcement.fit) OR autoplot(qcement.fit) isnt compatible with garchfit or doesn't seem to run in Rmarkdown at all.

** Results**

The ETS and ARIMA models are somewhat similar with little differences. I am using a quarterly dataset rather than a daily dataset here, so i tried to create quarterly lags (h=4) to reflect the dataset but I think I might've done it incorreclty. The residual time series appears to be very similar in all models. We also see some similarities between the ACF graphs in each model.