SS4861 Assignment 4&5 Han Zhang

14-4-3 SS4861 Assignment 4&5! Han Zhang hzhan463@uwo.ca The subject I choose to model and forecast is Canada Household final consumption expenditure(Food and non-alcoholic beverages),which represented as HFCE below.! First of all ,I download the data table from CANSIM and choose the unadjusted one.After this ,I plot the this time series using R and I get figure 1 below.! From figure 1,it is easy to see that this time series in seasonal and non-stationary,thus we consider the differenced time series wt=zt-zt-1 and we get the figure 2 below,which confirms the seasonal pattern and the lack of constant variance.! �1 Household final consumption expenditure(Food and non-alcoholic beverages), quarterly Time z 1980 1985 1990 1995 2000 2005 2010 2015 10 00 0 15 00 0 20 00 0 25 00 0 Differenced HFCE, quarterly Time w 1980 1985 1990 1995 2000 2005 2010 2015 -2 00 0 -1 00 0 0 10 00 20 00 14-4-3 According to figure 2,we can find that the variation is not constant and it seems to increase with the level of the series.This type of heteroscedasticity (non-constant variance) can usually be improved by using a logarithmic transformation. In statistics, we usually use natural logarithms for such transformations. ! Figure 3 below show the plot of wt=log(zt)-log(zt-1).! From the figure 3 above we can find that the variance is approximately constant.! Figure 4 below show the SACF of wt=log(zt)-log(zt-1).For seasonal time series with period s, the lags shown as 1, 2, 3, ... correspond to k =s, 2s, 3s, ... So our lag k corresponds on the plot to k/s in general.We look at the seasonal lags k = 0, s, 2 s, ... to select the seasonal component. In Figure 4 these seasonal lags correspond to 0, 1, 2, ... . The plot clearly shows that the SACF does not damp out quickly at the seasonal lags, hence we need ds>=1.! ! ! ! ! ! ! �2 Differenced log(HFCE), quarterly Time w 1980 1985 1990 1995 2000 2005 2010 2015 -0 .1 0 -0 .0 5 0. 00 0. 05 0. 10 14-4-3 ! ! �3 0 1 2 3 4 5 -0 .5 0. 0 0. 5 1. 0 Lag A C F logged, d=1 0 1 2 3 4 5 -0 .2 0. 0 0. 2 0. 4 0. 6 0. 8 1. 0 Lag A C F logged,d=1,ds=1,s=4 14-4-3 ! In Figure 6 above the first two lags are show at positions on the horizontal axis corresponding to 0.25 and 0.5 and these lags correspond to very large values of the PACF.! After fitting the model, the output is summarized below.! Call:! arima(x = log(z), order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))! Coefficients:! ar1 ar2 ma1 sma1! -0.1676 -0.1683 -0.0374 -0.7667! s.e. 0.4612 0.1199 0.4668 0.0609! sigma^2 estimated as 0.0001401: log likelihood = 381.46, aic = -752.92! The parameters are all statistically significant. Instead our goal should always be to fit the simplest possible model that passes the portmanteau diagnostic check.! The model diagnostic checks are shown in Figure 7 below and there is no evidence of lack of fit.! ! �4 1 2 3 4 5 -0 .2 0. 0 0. 2 0. 4 0. 6 Lag P ar tia l A C F logged, d=1, ds=1, s=4 14-4-3 Then we make forecasts and their standard errors were computed using the R predict() function. ! The results are displayed below:! > out2! $pred! Qtr1 Qtr2 Qtr3 Qtr4! 2014 10.03364 10.10563 10.11736 10.14195! 2015 10.06662 10.13977 10.15131 10.17574! ! $se! Qtr1 Qtr2 Qtr3 Qtr4! 2014 0.01183552 0.01511978 0.01702392 0.01902693! 2015 0.02218750 0.02463615 0.02671621 0.02871150! The above results are in the log domain. We need to transform back to the original data domain to get the forecasts for the original series. ! ! �5 Standardized Residuals Time 1980 1985 1990 1995 2000 2005 2010 2015 -3 -1 1 3 0 1 2 3 4 5 -0 .2 0. 2 0. 6 1. 0 Lag A C F ACF of Residuals 2 4 6 8 10 0. 0 0. 4 0. 8 p values for Ljung-Box statistic lag p va lu e 14-4-3 > zf! Qtr1 Qtr2 Qtr3 Qtr4! 2014 22780.09 24480.50 24769.31 25386.05! 2015 23543.74 25330.57 25624.53 26258.26! > zLo! Qtr1 Qtr2 Qtr3 Qtr4! 2014 22257.73 23765.67 23956.47 24456.76! 2015 22541.83 24136.50 24317.26 24821.40! > zUp! Qtr1 Qtr2 Qtr3 Qtr4! 2014 23314.71 25216.83 25609.73 26350.64! 2015 24590.19 26583.72 27002.08 27778.30! The resulting forecasts and their lower and upper 95% prediction limits are shown below in Figure 8 ,which contains the forecasts and their 50% prediction intervals are shown. ! �6 Time m ill io ns d ol la rs 2012 2013 2014 2015 22 00 0 23 00 0 24 00 0 25 00 0 26 00 0 27 00 0 28 00 0 Canada Household final consumption expenditure, observed and forecast 95% prediction interval shown in shaded region 14-4-3 Appendix R script! > setwd("/Users/zhanghan/Downloads")! > z<-ts(scan("ii.csv",sep=","),start=c(1978,1),freq=4)! Read 144 items! > plot(z, main="Household final consumption expenditure(Food and non-alcoholic beverages), quarterly")! > w<-diff(z)! > plot(w, main="Differenced HFCE, quarterly")! > #HFEC represents Household final consumption expenditure! > w<-diff(log(z))! > plot(w, main="Differenced log(HFCE), quarterly")! > acf(w, main="logged, d=1")! > w<-diff(z,lag=4,differences=1)! > acf(w,main="logged,d=1,ds=1,s=4")! > acf(w, main="logged, d=1, ds=1, s=4",type="partial")! > #ok try model! > out1 <- arima(log(z), order=c(2,1,1), seasonal=list(order=c(0,1,1),period=4))! > out1! Call:! arima(x = log(z), order = c(2, 1, 1), seasonal = list(order = c(0, 1, 1), period = 4))! Coefficients:! ar1 ar2 ma1 sma1! -0.1676 -0.1683 -0.0374 -0.7667! s.e. 0.4612 0.1199 0.4668 0.0609! ! sigma^2 estimated as 0.0001401: log likelihood = 381.46, aic = -752.92! > tsdiag(out1)! > #prediction! > out2 <- predict(out1, n.ahead=8)! > out2! �7 14-4-3 $pred! Qtr1 Qtr2 Qtr3 Qtr4! 2014 10.03364 10.10563 10.11736 10.14195! 2015 10.06662 10.13977 10.15131 10.17574! $se! Qtr1 Qtr2 Qtr3 Qtr4! 2014 0.01183552 0.01511978 0.01702392 0.01902693! 2015 0.02218750 0.02463615 0.02671621 0.02871150! > #forecasts and their prediction interval in log domain! > zfRaw <- out2$pred! > zfRawUp <- zfRaw + 1.96*out2$se! > zfRawLo <- zfRaw - 1.96*out2$se! > #need inverse log transformation! > zf <- exp(zfRaw)! > zUp <- exp(zfRawUp)! > zLo <- exp(zfRawLo)! > zf! Qtr1 Qtr2 Qtr3 Qtr4! 2014 22780.09 24480.50 24769.31 25386.05! 2015 23543.74 25330.57 25624.53 26258.26! > zLo! Qtr1 Qtr2 Qtr3 Qtr4! 2014 22257.73 23765.67 23956.47 24456.76! 2015 22541.83 24136.50 24317.26 24821.40! > zUp! Qtr1 Qtr2 Qtr3 Qtr4! 2014 23314.71 25216.83 25609.73 26350.64! 2015 24590.19 26583.72 27002.08 27778.30! > #to concatenate time series in R, it is easiest just to convert them! �8 14-4-3 > # to a vector, concatentate the vectors, and convert back to ts object! > z2 <- c(rev(rev(as.vector(z))[1:8]))! > zLast2 <- ts(c(z2, as.vector(zf)), start=2012, freq=4)! > #since we want to show 95% prediction limit, we need to use the ylim argument in plot! > #choose the factors 0.97 and 1.02 by the trial-and-error method! > ylim <- c(0.97*min(zf, zLo), 1.005*max(zf,zUp))! > plot(zLast2, ylim=ylim,ylab="millions dollars")! > points(zf, col="red", cex=1.5, pch=16)! > points(ts(z2,start=2012,freq=4), col="black", cex=1.5, pch=16)! > legend(x=2012.25, y=56*10^4, legend=c("observed","forecast"),! + pch=c(16,16),col=c("black","red") )! > #lets show the confidence region as colored region! > #we use rgb function with low alpha setting in polygon()! > x <- as.vector(time(zLo))! > ylo <- as.vector(zLo)! > yup <- as.vector(zUp)! > xx <- c(x, rev(x))! > yy <- c(ylo,rev(yup))! > polygon(x=xx, y=yy, border=FALSE, col=rgb(0.5,0.5,0,0.4))! > title(main="Canada Household final consumption expenditure, observed and forecast\n95% prediction interval shown in shaded region") �9