A. Get to know your data set.



    Variable Name:

    Bicol Pork Lean Meat

    Define your time series with appropriate unit of measurement:

    Monthly retail price of pork lean meat in Bicol Region from January 1995 to July 2018 where:

    RP is used to denote the retail price (in Php/kg) of Pork Lean Meat in Bicol region; and

    Time is used to denote the monthly period. For model building, the data from January 1995 to December 2016 is used. While for forecast validation, the data from January 2017 to July 2018 is used.

    Sample size:

    283 observations (the first 264 observations for model building; the remaining 19 observations for forecast evaluation)

    Source:

    Philippine Statistics Authority (PSA)




B. Characterizing the time series

    1. Graph your time series and through ocular inspection, identify what components of the time series are present, i.e., trend, drift, cycle or seasonality.

      Solution:
      As a practice in time series, the original data sets are divided into two sets - the first one (larger set) is for model building and the second set (validation set) is for forecast evaluation. In this problem, the Bicol pork lean meat retail price data covering January 1995 to December 2016 is used in the model building.

      Using the R code below, we can obtain the graph of the given time series data as shown in Figure 1.

      library(astsa)
library(readxl)
      #Load the data file
DataLoc = "D:/MSAM/Math 215 - Time Series Analysis/Bangi_FinalExam - Data.xlsx"
BPLM_Data <- read_excel(DataLoc, sheet="BICOL-PORK", range="A1:D265")

#Plot the time series data
Data = ts(BPLM_Data$RP, frequency = 12, start = c(1995,1))
plot.ts(Data, main="Bicol Pork Lean Meat Retail Price", 
        ylab="RP", xlab="Time", col=4)

      Figure 1. Bicol pork lean meat retail prices from January 1995 to December 2016


      Through ocular inspection, we can see an upward trend in the retail price of pork lean meat in the Bicol region from January 1995 to December 2016. When we look at it closely, we can observe that there is a presence of multiple trends and drifts. This is specified in Figure 2.

      Also, if we look at the overall data in Figure 3, we can somehow see a cycle. Note that a cycle is a periodicity in the given time series data with interval more than annual. Details of which is presented in part (b).


      Figure 2. Presence of trends and drifts in the Bicol pork lean meat retail prices data



      Figure 3. Presence of periodicity in the Bicol pork lean meat retail prices data


    1. Explain one by one why you think these components are present in your time series.

      Solution:

      Trend
      As shown in Figure 2, there are three trends identified where Trend 1 and Trend 3 are somehow similar since their trend lines are showing almost a parallel increase in the retail prices of Bicol pork lean meat. Trend 2 is showing a period of price increase followed by a period of price decrease. Moreover, the rate of increase of the trends are not the same. Trend 2 appears to be steeper compared to the other trends.


      Drift
      The change in the rate of the increasing trends is due to the sudden change in the mean level of the time series somewhere around 2001 and 2011. These suggest that drifts are present in the time series data.


      Periodicity
      If we focus on the data inside the red box in Figure 3, we can see that in certain period of months, the price of pork lean meat in the Bicol region increase then it drops. This appears to have two and half cycles of such price behavior. Moreover, as mentioned above, Trend 1 and Trend 3 are somehow parallel (refer to the encircled data). The behavior of the data covering these trends also appear to be near identical.




C. Stationarizing the time series

    Definition: Given \(\small Y_t\)is a time series, it is weakly stationary if its mean and variance are constant with respect to time. There are two ways by which the trend and periodicity can be removed from the series. These methods are:

    1. Detrending
      Estimate the trend component of the series by fitting the linear regression model as follows: \(\small Y_t=\alpha + \beta_1 time + \mu_t\)

      1. Estimate the trend line and superimpose it with the original data \(\small Y_t\).

        Solution:
        Note that from the given linear regression model \(\small Y_t=\alpha + \beta_1 time + \mu_t\), we let \(\small \alpha\) be the intercept (the value of the trend at time \(\small t=0\)) and \(\small \beta_1\) be the slope.

        Using the R code below, we can obtain the fit summary.

        #Original data with trend line
trend = time(Data)
summary(fit <- lm(Data~trend, na.action=NULL))
        ## 
## Call:
## lm(formula = Data ~ trend, na.action = NULL)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.6728  -5.2315   0.0966   4.5001  17.7681 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.145e+04  1.305e+02  -87.75   <2e-16 ***
## trend        5.780e+00  6.508e-02   88.81   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.715 on 262 degrees of freedom
## Multiple R-squared:  0.9679, Adjusted R-squared:  0.9677 
## F-statistic:  7888 on 1 and 262 DF,  p-value: < 2.2e-16
        coef(fit)
        ##   (Intercept)         trend 
## -11454.517855      5.779787

        Looking at the result, we have obtained the intercept \(\small \alpha\) and slope \(\small \beta_1\) values equal to -11454.517855 and 5.779787, respectively. Thus, the estimated trend can be written as \[\small \hat{\mu}_t=-11454.517855+5.779787t\] where \(\small t\) is time expressed as monthly units. Also, we can observe that the slope \(\small \beta_1\) is positive which implies that the trend is increasing. The larger the absolute value of \(\small \beta_1\), the steeper the trend’s slope becomes.

        Plotting the graph of the estimated trend line and superimpose it with the original data \(\small Y_t\), as shown in Figure 4, we can use the R code below.

        plot.ts(Data, ylab="RP", lwd=2, col=4,
        main="Linear Trend of the Bicol Pork Lean Meat Retail Price Data")
abline(fit, col=2)

        Figure 4. Linear trend of the Bicol pork lean meat retail prices data


      1. Estimate the detrended series (residual component \(\small \mu_t\)).

        Solution:
        Let \(\small x_t\) be the observations at time \(\small t\). To obtain the detrended series, we subtract \(\small \hat{\mu}_t\) obtained in part (a) from the observations. Hence, the estimated detrended series is \[\small \hat{y}_t=x_t+11454.517855-5.779787t\] Using the R code below, we obtain the graph shown in Figure 5 showing the detrended series.

        plot.ts(resid(fit), main="Detrended time series data", 
        ylab="Residuals", col=8)

        Figure 5. Detrended time series data


      Is the error component a white noise (i.e., stationary with respect to the mean and variance of the series)? Explain your answer.


      Solution:
      As shown in Figure 5, the series appears to be stationary with respect to the mean of zero. However, if we check the autocorrelation function (ACF) of the detrended series shown in Figure 6, we can observe that it slowly decays. Note that ACF describes how well the present value of the series is related with its past values. The detrended ACF has shown lags 0 to 14, lags 25 to 30, lags 45 to 51, lags 62 to 99, lags 119 to 130, lags 152 to 161, and lags 189 to 190 to be outside of the horizontal dashed lines, which are placed at zero plus and minus two approximate standard errors of the sample autocorrelations, namely \(\small \pm2 / \sqrt{n}\). This is not what we expect from a white noise process.

      Below is the R code to generate the ACF plot of the detrended series.

      acfpl_detrended <- acf(resid(fit), 263, plot=FALSE)
acfpl_detrended$lag <- acfpl_detrended$lag * 12
plot(acfpl_detrended, xlab="Lag (months)",
     main="ACF of Detrended Series", col=2)

      Figure 6. Sample ACF of the detrended series


    1. Differencing
      Difference the series to stationarize the original series. Obtain the first difference of the series by using the following formula: \(\small \Delta Y_t=Y_t-Y_{t-1}\)


      1. Graph the 1st differenced series.

        Solution:
        Using the R code below, we obtain the graph shown in Figure 7.

        plot.ts(diff(Data), main="First Difference", ylab="Difference", col=8)

        Figure 7. First difference of Bicol pork lean meat retail prices





      1. Was the series stationarized with respect to the mean and variance after first differencing? If yes, explain your answer. If no, do you think the second difference will be needed? Explain your answer.

        Solution:
        As shown in Figure 7, the first differenced series appears to be stationary with respect to the mean of zero (mean is constant over the series). It is also apparent that the first differencing removed the trend in the series.

        If we plot the ACF of the first differenced series as shown in Figure 8, we can further observe that it decays to zero faster.

        The R code to generate the first difference ACF is found below.

        acfpl_diff <- acf(diff(Data), 263, plot=FALSE)
acfpl_diff$lag <- acfpl_diff$lag * 12
plot(acfpl_diff, xlab="Lag (months)",
     main="ACF of First Difference", col=2)

        Figure 8. First differencing ACF of Bicol pork lean meat retail prices


    1. Which procedure (i.e., detrending or differencing) is more appropriate in stationarizing your time series. Justify your answer.

      Solution:
      A stationarized time series data is characterized by not having any cyclical or seasonal patterns, nor any decreasing or increasing trend, and randomly fluctuates around the mean. Detrending is a procedure that removes any drift in the series but does not remove its increasing/decreasing variance. On the other hand, differencing helps to stabilize the mean.

      Comparing the results in C.1 and C.2, it is evident that the differencing procedure is more appropriate in stationarizing the time series data. Detrending the data still exhibits non-stationarity with its slow decay; while getting its first difference resulted to a stationary series.

      This observation can be verified using the Augmented Dickey-Fuller test in part D.



D. Testing for stationarity using Augmented Dickey-Fuller test

    1. Execute the following tests for the level form of the series:

      1. Unit root test with a drift

        Solution:
        We establish the hypothesis for the level form of the series using the Augmented Dickey-Fuller (ADF) Test with drift as follows:

          \(\small H_0\): The Bicol lean pork retail price is non-stationary at level form.
          \(\small H_a\): The Bicol lean pork retail price is stationary at level form.


        Using the ‘adfTest’ function with type set to “c” in the ‘fUnitRoots’ library of R, we obtain the following table:

        Table 1. Test for stationarity of the level form (Unit root test with drift)


        As shown in Table 1, the p-value of ADF Test for the level form data (0.7569) is higher than 0.05. Therefore, we failed to reject \(\small H_0\). This implies that the Bicol lean pork retail price is non-stationary at level form when unit root test with drift is used.

        Below is the R code used to generate the above result.

        library(fUnitRoots)
        adfTest(Data, type = "c")


      1. Unit root test with trend

        Solution:
        We establish the hypothesis for the level form of the series using the Augmented Dickey-Fuller (ADF) Test with trend as follows:

          \(\small H_0\): The Bicol lean pork retail price is non-stationary at level form.
          \(\small H_a\): The Bicol lean pork retail price is stationary at level form.


        Using the ‘adfTest’ function with type set to “ct” in the ‘fUnitRoots’ library of R, we obtain the following table:

        Table 2. Test for stationarity of the level form (Unit root test with trend)


        As shown in Table 2, the p-value of ADF Test for the level form data (0.1831) is slightly higher than 0.05. Therefore, we failed to reject \(\small H_0\). This implies that the Bicol lean pork retail price is non-stationary at level form when unit root test with trend is used.

        Below is the R code used to generate the above result.

        adfTest(Data, type = "ct")


      Discuss the result of each test, i.e., is the series stationary or not. Defend your answer. Which is the better model of the stationary series?


      Solution:
      From D.1.i and D.1.ii, both results indicate that the time series at level form is not stationary. If we compare the results between these unit root tests, it is apparent that the unit root test with trend is a better one.



    1. Execute the following tests for the first difference of the series:

      1. Unit root test with a drift: \(\small \Delta Y_t=\alpha+\rho \Delta Y_{t-1}+\epsilon_t\)

        Solution:
        We establish the hypothesis for the first difference of the series using the Augmented Dickey-Fuller (ADF) Test with drift as follows:

          \(\small H_0\): The Bicol lean pork retail price is non-stationary at the first difference of the series.
          \(\small H_a\): The Bicol lean pork retail price is stationary at the first difference of the series.


        Using the ‘adfTest’ function with type set to “c” in the ‘fUnitRoots’ library of R, we obtain the following table:

        Table 3. Test for stationarity of the transformed series (Unit root test with drift)


        As shown in Table 3, the p-value of ADF Test for the first difference is lower than 0.05. Therefore, we reject \(\small H_0\) and accept \(\small H_a\). This implies that the Bicol lean pork retail price is stationary at the first difference of the series when unit root test with drift is used.

        Below is the R code used to generate the above result.

        adfTest(diff(Data), type = "c")


      1. Unit root test with trend: \(\small \Delta Y_t=\alpha+\rho \Delta Y_{t-1}+\beta time+\epsilon_t\)

        Solution:
        We establish the hypothesis for the first difference of the series using the Augmented Dickey-Fuller (ADF) Test with trend as follows:

          \(\small H_0\): The Bicol lean pork retail price is non-stationary at the first difference of the series.
          \(\small H_a\): The Bicol lean pork retail price is stationary at the first difference of the series.


        Using the ‘adfTest’ function with type set to “c” in the ‘fUnitRoots’ library of R, we obtain the following table:

        Table 4. Test for stationarity of the transformed series (Unit root test with trend)


        As shown in Table 4, the p-value of ADF Test for the first difference is lower than 0.05. Therefore, we reject \(\small H_0\) and accept \(\small H_a\). This implies that the Bicol lean pork retail price is stationary at the first difference of the series when unit root test with trend is used.

        Below is the R code used to generate the above result.

        adfTest(diff(Data), type = "ct")


      Discuss the result of each test, i.e., is the series stationary or not. Defend your answer. Which is the better model of the stationary series?


      Solution:
      From D.2.i and D.2.ii, we see that both the results at the first difference of the time series data are stationary having p-value less than 0.01. Note further that since it is already assumed that the first differencing eliminates the trend of any given series, it is most likely that the difference series would not need the type set to “ct”.



E. ARIMA Prediction

    1. Obtain the ACF of the time series at level form. Does the ACF plot suggest stationarity? Compare with your result in D.2.

      Solution:
      Looking at the ACF of the time series at level form, as shown in Figure 9, the monthly retail price of Bicol pork lean meat slowly decays and does not cut off at fewer lags. These suggest non-stationarity. When we further performed ADF Test in D.1, we had verified that the generated p-value is higher than 0.05 which implies that the time series data at level form is indeed non-stationary. In order to stationarize this, we performed the first differencing in C.2 and verified using ADF Test in D.2.

      Below is the R code to generate the ACF of the original data.

      acfpl_orig <- acf(Data, 263, plot=FALSE)
acfpl_orig$lag <- acfpl_orig$lag * 12
plot(acfpl_orig, xlab="Lag (months)",
     main="ACF of Bicol Pork Lean Meat Retail Price", col=2)

      Figure 9. ACF plot of the original data


    1. Determine the tentative models. Discuss how you were able to identify these models. Compare their AIC values and Discuss which of these models is the final model.

      Solution:
      Using the R codes below, we can generate both the ACF and PACF plots of the first difference of the series.

      acfpl_diff <- acf(diff(Data), plot=FALSE)
acfpl_diff$lag <- acfpl_diff$lag * 12
pacfpl_diff <- pacf(diff(Data), plot=FALSE)
pacfpl_diff$lag <- pacfpl_diff$lag * 12

par(mfrow = c(1,2))
plot(acfpl_diff, xlab="Lag (months)",
     main="ACF of the First Difference", col=2)
plot(pacfpl_diff, xlab="Lag (months)",
     main="PACF of the First Difference", col=2)
      Figure 10. Sample ACF (left) and Partial ACF (right) of the first difference series


      As observed in Figure 10, we note the following:

        • The ACF cuts off at lag 2 while the PACF appears to be a mixture of exponential decay and damped sinusoid. With this observation, we can consider the model ARIMA (0,1,2).

        • The ACF appears to be a mixture of exponential decay and damped sinusoid while the PACF cuts off at lag 1. With this observation, we can consider the model ARIMA (1,1,0).

        • Both the ACF and PACF tails offs suggesting a model with both AR and MA terms. With this observation, we can consider the models ARIMA (1,1,1) and ARIMA (1,1,2).

      Table 5 shows the tentative ARIMA models to forecast the monthly pork lean meat retail prices in Bicol region together with its respective Akaike’s Information Criterion (AIC) values. Note that the AIC is used to estimate the model quality. Among the tentative models, ARIMA (0,1,2) is selected since it has the least AIC value of 1042.824.

      Table 5. Tentative ARIMA Models


      Below are the R codes used to generate the AIC values of the tentative models.

      library(forecast)
      #Tentative Models
M1 = Arima(Data, order=c(0,1,2))
M2 = Arima(Data, order=c(1,1,0))
M3 = Arima(Data, order=c(1,1,1))
M4 = Arima(Data, order=c(1,1,2))
      #Get the AIC
M1$aic
M2$aic
M3$aic
M4$aic


    1. Fit the ARIMA model by estimating the model parameters. Are the estimates of the parameters significant?

      Solution:
      The estimates of the model’s coefficients for ARIMA (0,1,2) are presented in Table 6 along with the standard error, z-values, and p-values. As observed, the MA (1) has an estimate value of 0.280336 while MA (2) has an estimate value of 0.282498. Both model terms have p-values less than the significance level of 0.05. Hence, both coefficients are significantly different from zero.

      Table 6. Model Estimates of ARIMA (0,1,2)


      Below are the codes used to generate the model estimates.

      library(lmtest)
      #Set the best model
Model = M1
      coeftest(Model)


    1. Perform residual analysis. Generate residual plots (residual vs time, residual vs fitted, normal qq plot, acf and pacf plots). Discuss the results thoroughly.

      Solution:
      Residual analysis is performed to check the adequacy of the model. First, we plot the residuals against time using the R code below.

      #Residual Analysis
res = Model$residuals
fit = Model$fitted
      #Plot of residuals vs. time
plot.ts(res, type="p", ylab="Residuals", main="Residuals vs. Time", col=4)
abline(0,0)
      Figure 11. Plot of Residuals vs. Time


      Figure 11 shows the plot of residuals vs. time. The residuals do not show any strong patterns and fluctuate randomly around zero. This suggests that the residuals are independent against time.

      Next, we plot the residuals against the fitted values. This can be generated using the R code as follows:

      #Plot of residuals vs. fitted values
plot.ts(as.numeric(fit), as.numeric(res), xlab="Fitted Value", 
     ylab="Residuals", main="Residuals vs. Fitted Values", col=4)
abline(0,0)
      Figure 12. Plot of Residuals vs. Fitted Values


      Figure 12 shows the plot of residuals vs. fitted values. It can be observed that the plot has no visible pattern. Hence, we can assume that the constancy of error variance is met.

      To generate the normal Q-Q plot of residuals, we can use the R code below.

      #Normal Q-Q plot of residuals
qqnorm(res, main="Normal Q-Q Plot for Residuals", col=4)
qqline(res, col=1)
      Figure 13. Normal Q-Q Plot of Residuals


      We can observe in Figure 13 that some residuals seem to deviate from the theoretical line of normality. This suggests that the residuals are not normally distributed. This can be further verified using the Shapiro test for normality which generates a p-value less than 0.05 (see the result below).

      #Check for Normality of Residuals
shapiro.test(res)
      ## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.95165, p-value = 1.127e-07

      Plotting the ACF and PACF plots of residuals, we use the following R code:

      #ACF and PACF of residuals
acfpl_res <- acf(res, plot=FALSE, main="ACF of Residuals", col=2)
acfpl_res$lag <- acfpl_res$lag * 12
pacfpl_res <- pacf(res, plot=FALSE, main="PACF of Residuals", col=2)
pacfpl_res$lag <- pacfpl_res$lag * 12
par(mfrow = c(1,2))
plot(acfpl_res, xlab="Lag (months)", main="ACF of Residuals", col=2)
plot(pacfpl_res, xlab="Lag (months)", main="PACF of Residuals", col=2)
      Figure 14. Sample ACF (left) and Partial ACF (right) of Residuals


      In Figure 14, it can be observed that most of the ACFs and PACFs are within the acceptable limits although there are spikes in both ACF and PACF at lag 4. Thus, the residuals can still be considered a white noise process.

      We can further perform the Ljung-Box Test to formally test the lack of fit of the model. That is, we define the following:

        \(\small H_0\): The ARIMA (0,1,2) model does not show lack of fit.
        \(\small H_a\): The ARIMA (0,1,2) model does show a lack of fit.


      Table 7 shows the Ljung-Box test statistics value of 33.6 with a p-value of 0.05391 which is slightly greater than the significance level of 0.05. This implies that the null hypothesis cannot be rejected. Therefore, we can conclude that the ARIMA (0,1,2) model fits the data.

      Table 7. Ljung-Box Test for lack of fit of ARIMA (0,1,2)


      Below is the R code used to obtain the Ljung-Box Test result.

      checkresiduals(Model)


    1. Generate the one-step ahead forecast. Plot the ACF and PACF plots of forecast errors and test its normality. Discuss thoroughly your observations.

      Solution:
      Table 8 summarizes the actual retail prices, the generated one-step ahead forecasted values, and the corresponding forecast errors of ARIMA (0,1,2) model. It can be observed that the number of positive forecast errors is larger than the number of negative forecast errors. This implies that most of the actual retail prices are larger than the forecasted values. Moreover, it is interesting to note that the forecast errors do not exceed \(\small \pm 7\) (in Php/kg).

      Table 8. One-Step Ahead Forecast and Forecast Errors (in Php/kg)


      This is obtained using the codes below.

      #Validation - One-Step Ahead Forecast
BPLM_DataVal <- read_excel(DataLoc, sheet="BICOL-PORK-VAL")
DataVal = ts(BPLM_DataVal$RP)
NewModel = Arima(c(Data, DataVal), model=Model)
Onestep = fitted(NewModel)[265:283]
Onestep
#Forecast Errors
FError = DataVal - Onestep
FError


      Moreover, Figure 15 shows the graphical result of the one-step ahead forecast values (red) against the actual (blue) retail prices of Bicol pork lean meat from January 2017 to July 2018. It can be generated using the R code below.

      plot.ts(DataVal, main="Actual vs. Forecasted",
        ylab="RP", lwd=2, col=4, xaxt="n")
axis(1, at=1:19, labels=BPLM_DataVal$Year)
lines(Onestep, col=2)
      Figure 15. Actual vs. Forecasted Retail Price of Bicol Pork Lean Meat


      On the other hand, the ACF and PACF plots of the forecast errors are presented in Figures 16. Based on the plots, most of the ACFs and PACFs are within the acceptable limits indicating that the values are not significantly different from zero. Thus, the forecast errors generally behave like a white noise.

      The plots are generated using the R code below.

      par(mfrow = c(1,2))
acf(FError, 18, main="ACF of Forecast Error", col=2)
pacf(FError, 18, main="PACF of Forecast Error", col=2)
      Figure 16. Sample ACF and Partial ACF of Forecast Errors


      To formally check if the forecast errors are normally distributed, the Shapiro-Wilk Test is applied. Below is the R code to perform the test.

      shapiro.test(FError)
      Table 9 shows the test statistics value and its p-value.

      Table 9. Shapiro Test for Normality of Forecast Errors


      Since the p-value (0.5858) is greater than significance level of 0.05, the null hypothesis that the error terms are normally distributed cannot be rejected. Hence, the forecast errors can be considered as a Gaussian white noise.



    1. Generate a forecast for the next ten periods.

      Solution:
      Since the diagnostics of the model and the forecast evaluation suggest that the ARIMA (0,1,2) model is adequate for forecasting, we can now generate forecast values for the next ten months. The forecasted retail prices of Bicol pork lean meat covering the period from August 2018 to May 2019 are shown in Table 10 together with its lower and upper 95% confidence interval. Please note that even if the forecasted RP appears to be the same from September 2018 to May 2019, we note that its lower and upper 95% confidence intervals differ. This implies that actual retail prices may fall within these thresholds. Once a new actual data is added, the forecasted RP will then be updated.

      Table 10. Forecasted Retail Price of Pork Lean Meat in Bicol Region (in Php/kg)


      These values are generated using the R code below.

      Forecasts <- forecast(c(Data, DataVal), model=Model, h=10)
fc = as.data.frame(Forecasts)
print(fc)

      Figure 17 further shows the plot of the original data together with the forecasted values and confidence intervals for the next ten periods.

      xlabels=c(BPLM_Data$Year,BPLM_DataVal$Year,NA,NA,NA,NA,NA,2019,NA,NA,NA,NA)
plot(Forecasts, col=4, fcol=2, ylab="RP", xlab="Month", xaxt="n")
axis(1, at=1:293, labels=xlabels)

      Figure 17. Plot of Original Data with the Forecasted Values and Confidence Intervals



— Nothing follows —