Forecasting Principles Practices

Assignment

Forecasting Principles and Practices 

Chapter 7 Exercises

Data set books contains the daily sales of paperback and hardcover books at the same store.  The task is to forecast the next four days’ sales for paperback and hardcover books (data set books).

7.2          (a) Apply Holt’s linear method to the paperback and hardback book series and compute four-day            forecasts in each case.

(b) Compare the SSE measures of Holt’s method for the two series of those simple exponential                               smoothing in the previous question. Discuss the merits of the two forecasting methods for those            two data sets.

(c) Compare the forecasts for the two series using both methods.  Which do you think is best?

(d) Calculate a 95% prediction interval for the first forecast for each series using both methods,                 assuming normal errors.  Compare your forecasts with those produced by R.

7.3 For this exercise, use the price of a dozen eggs in the United States from 1900 – 1993 (data set eggs).  Experiment with the various options in the holt () function to see how much the forecasts change with damped or exponential trend.  Also try change the parameter values for α and to see how they affect the forecasts.  Try to develop an intuition of what each parameter and the argument is doing to the forecasts.

[Hint:  use h=100 when calling holt () so you can clearly see the differences between the various options when plotting the forecasts.]

(b) Which model gives the best RMSE? 

Chapter 8 Exercises

8.5 Use R to simulate and plot some data from simple ARIMA models.

(a) Generate data from an AR(1) model with  = 0.6 and  =1. Start with  = 0.

(b) Produce a time plot for the series. How does the plot change as you change   ?

(c) Generate data from an MA (1) model with  = 0.6 and  =1. Start with = 0.

(d) Produce a time plot for the series.  How does the plot change as you change  ?

(e) Generate data from an ARIMA (1, 1) model with  = 0.6 and  = 0.6 and  =1.  Start with  = 0 and = 0.

(f) Generate data from an AR (2) model with  = – 0.8 and   = 0.3 and  =1. Start with   =                  =  0.

(g) Graph the latter two series and compare them.

8.6 Consider the number of women murdered each year (per 100,000 standard population) in the United States (datasetwmurders).

(a)  By studying appropriate graphs of the series in R, find and appropriate ARIMA (p, d, q) model            for these data.

(b) Should you include a constant in the model? Explain.

(c) Write this model in terms of a backshift operator.

(d) Fit the model using R and examine the residuals.  Is the model satisfactory?

(e) Forecast three times ahead.  Check your forecasts by hand to make sure you know how they             have been calculated.

(f) Create a plot of the series with forecasts and prediction intervals for the next three periods                 shown.

(g) Does auto.arima give the same model you have chosen?  In not, which model do you think is              better?

(h) Find the latest data from:

https://catalog.data.gov/dataset?q=Murders&sort=score+desc%2C+name+asc&as_sfid=AAAAAAU2omFQEqfdZ_e56PcYRNBl0Jfl26t4g56JZTrb-a8NSOgssWqvXWSk0lVPO4YrZRJID4N-DaojOmM_BNHhTY7O2VcagygjnfGC4ddASpswBduRIP0Dilrjah1gmwWLfo0%3D&as_fid=d7b23e6858d59ffba454b1e14e0f7a790f56e9cd&ext_location=&ext_bbox=&ext_prev_extent=-142.03125%2C-3.162455530237848%2C-59.0625%2C55.57834467218206

8.7 Consider the quarterly number of international tourists to Australis for the period 1999-2010 (data set austourists).

(a) Describe your time plot.

(b) What can you learn from the ACF graph?

(c) What can you learn from your PACF graph?

(d) Produce plots of the seasonally differenced data  .  What model do these graphs               suggest                ?

(e)  Does auto.arima give the same model that you chose?  If not, which model do you think is better?

(f) Write the model in terms of the backshift operator, and then without using the backshift       operator.

8.8 Consider the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period 1985 – 1996). (Data set usmelec). In general there are two peaks per year:  in mid-summer and mid-winter.

(a) Examine the 12-month moving average of this series to see what kind of trend is involved.

(b) Does the data need transforming?  If so, find a suitable transformation.

(c) Are the data stationary?  If not, find an appropriate differencing which yields stationary data.

(d) Identify a couple of ARIMA models that might be useful in describing the time series.  Which               of your models is the best according to AIC values?

(e) Estimate the parameters of your best model and do diagnostic testing on the residuals.  Do                 the residuals resemble white noise? If not, try to find another ARIMA model which fits better.

(f) Forecast the next 15 years of generation of electricity by the U.S. Electric Industry.  Get the latest figures from https://catalog.data.gov/dataset?q=Electricity&sort=score+desc%2C+name+asc&as_sfid=AAAAAAXtASDDEhyS1_7pDTIJf-sWEHU66YEtjZdYG8QD_LDTKOM74IwGTnmUp86zk6DH_v7cS-x3ImrOOEgbwPTt7UFk-UuCW7QFrJqQDcHhZA_yXeyEYGX0XcqJ2QFjxOhDGZE%3D&as_fid=f0425c4fe86bd75deaed5fdb5c571629254ad9a0&ext_location=&ext_bbox=&ext_prev_extent=-142.03125%2C8.754794702435617%2C-59.0625%2C61.77312286453146

to check the accuracy of the forecasts.

(g) How many years of forecasts do you think are sufficiently accurate to be suitable?

Solution 

title: “timeseries16102017″

output: word_document

“`{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

“`

## Load library

“`{r, echo = FALSE}

library(fpp)

library(ggplot2)

library(ZIM)

library(quantmod)

“`

##Chapter 7 Exercises

Dowloads dataset

“`{r }

df<- books

class(df)

“`

#a

“`{r }

fit1 <- holt(df[,”Paperback”], alpha=0.8, beta=0.2, initial=”simple”,  h = 4)

fit1$model$states

fitted(fit1)

fit1$mean

plot(fit1, type = “o”, main = “Forecast Paperback”)

fit2 <- holt(df[,”Hardcover”],alpha=0.8, beta=0.2, initial=”simple”, h = 4)

plot(fit2, type = “o”,  main = “Forecast Hardcover”)

“`

#b

We built two models for forecasting Forecast Paperback and Hardcover. We see that the coefficient SSE is smaller in the second model, which predicts “Forecast Hardcover”. SSE for fit2 is 44383.94.

SSE is the sum of the squared differences between each observation and its group’s mean. It can be used as a measure of variation within a cluster. If all cases within a cluster are identical the SSE would then be equal to 0. If you look at the trend dynamics, then the trends of both variables are slightly different. It should be assumed that for the model Holt’s linear method the variable Paperback has a somewhat chaotic nature. Therefore, to reduce errors, it is necessary either to transform the data, or to try to apply other models.

“`{r }

autoplot(df)

fit1$model$SSE

fit2$model$SSE

“`

#c

I have slightly transformed the input parameters into models.I changed the alpha and beta ratios. It should be noted that SSE errors have decreased. Thus, through the adjustment of the model parameters, we can improve its accuracy. If you look at the forecast graph, you should pay attention that we see the continuation of the upward trend, but with a fairly wide confidence interval.

“`{r }

fit1.exp <-  holt(df[,”Paperback”], alpha=0.2, beta=0.8, initial=”simple”,  h = 4, exponential = T)

fit2.exp <- holt(df[,”Hardcover”],alpha=0.2, beta=0.8, initial=”simple”, h = 4, exponential = T)

fit1.exp$model$SSE

fit2.exp$model$SSE

plot(fit1.exp, type = “o”, main = “Forecast Paperback with exponential smoothing”)

plot(fit2.exp, type = “o”, main = “Forecast Hardcover with exponential smoothing”)

“`

#D

I calculated a 95% confidence interval for each model.I separately calculate the upper and lower confidence interval. If we analyze the trend, then this will be our levels of support and resistance. If you look at the dynamics of the time series, you can see that some models very clearly predict short-term “seasonal fluctuations,” while other models capture longer-term “seasonal fluctuations.” This means that setting up the model has a very important role in creating an accurate model.

Interestingly, the first non-modalized models show a more pronounced oscillation dynamics, although they have SSE large errors. Modified models have smaller SSE errors and it is visually evident that the forecast plot for these models is smoother.

“`{r }

fit.conf.h<- fitted(fit1) + 1.96*(sd( fitted(fit1)))

fit.conf.l<- fitted(fit1) – 1.96*(sd( fitted(fit1)))

autoplot(df[,”Paperback”], PI = F, series = “Paperback original”) +

autolayer(fitted(fit1), PI = F, series = “Holt`s linear”) +

autolayer(fit.conf.h, PI = F, series = “High confidence layer”) +

autolayer(fit.conf.l, PI = F, series = “Low confidence layer”)

fit2.conf.h <- fitted(fit2) + 1.96*(sd( fitted(fit2)))

fit2.conf.l <- fitted(fit2) – 1.96*(sd( fitted(fit2)))

autoplot(df[,”Hardcover”], PI = F, series = “Hardcover original”) +

autolayer(fitted(fit2), PI = F, series = “Holt`s linear”) +

autolayer(fit2.conf.h, PI = F, series = “High confidence layer”) +

autolayer(fit2.conf.l, PI = F, series = “Low confidence layer”)

fit1.exp.conf.h <- fitted(fit1.exp) + 1.96*(sd( fitted(fit1.exp)))

fit.exp.conf.l<- fitted(fit1.exp) – 1.96*(sd( fitted(fit1.exp)))

autoplot(df[,”Paperback”], PI = F, series = “Paperback original”) +

autolayer(fitted(fit1.exp), PI = F, series = “Holt`s linear”) +

autolayer(fit1.exp.conf.h, PI = F, series = “High confidence layer”) +

autolayer(fit.exp.conf.l, PI = F, series = “Low confidence layer”)

fit2.exp..conf.h <- fitted(fit2.exp) + 1.96*(sd( fitted(fit2.exp)))

fit2.exp.conf.l <- fitted(fit2.exp) – 1.96*(sd( fitted(fit2.exp)))

autoplot(df[,”Hardcover”], PI = F, series = “Hardcover original”) +

autolayer(fitted(fit2.exp), PI = F, series = “Holt`s linear”) +

autolayer(fit2.exp..conf.h, PI = F, series = “High confidence layer”) +

autolayer(fit2.exp.conf.l, PI = F, series = “Low confidence layer”)

“`

#7.3

On the chart, we see a clear downtrend with some volatility. I changed the parameters of the model through the adjustment of the alpha and beta coefficients and compared various metrics. The model with metrics aplpha = 0.9, beta = 0.1 showing good results.

“`{r }

df.eggs<- eggs

autoplot(eggs)

tunihg.model<- holt(eggs, alpha = 0.9, beta = 0.1, initial = ‘simple’, h = 100)

tunihg.model$model$SSE

round(accuracy(tunihg.model))

“`

##8.5

#a

We make the function for the simulation of the AR1 process

“`{r }

y <- ts(numeric(100))

e <- rnorm(100)

for(i in 2:100)

y[i] <- 0.6*y[i-1] + e[i]

“`

#b

Lags ACF and PACF fairly quickly approach zero.

We see a fairly volatile schedule. The graph ACF and ACF shows that the first lag is beyond the confidence interval. This dynamic is characteristic of the process AR1

“`{r }

autoplot(y)

tsdisplay(y)

“`

#c

We make the function for the simulation of the MA(1) process

“`{r }

y <- ts(numeric(100))

e <- rnorm(100)

sim.data.mal<- function(n.obs, theta, seed.nr){

set.seed(seed.nr)

y <- ts(numeric(n.obs))

e <- rnorm(n.obs)

for (i in 2:n.obs)

y[i] <- theta*y[i-1] + e[i]

return(y)

}

“`

#d

We see a fairly volatile trend dynamics. Lags have a cyclic shape, which is typical for the process MA. Also we see that the first   lags  in ACF ans PACF go beyond the confidence interval. This is a process MA(1)

“`{r }

plot(sim.data.mal(100, 0.6, 1), main=”time series”, xlab=”time”)

tsdisplay(sim.data.mal(100, 0.6, 1))

“`

#e

Now let’s make the function that will generate the process ARIMA(1,1).

We see a changeable trend. The first lags in ACF and PACF go beyond their confidence interval. In addition, there is a certaincyclicity in the graph.

“`{r }

sim.data.armal1 <- function(n.obs, phi, theta, seed.nr){

set.seed(seed.nr)

y <- ts(numeric(n.obs))

e <- rnorm(n.obs)

for (i in 2:n.obs)

y[i] <- phi*y[i-1] + theta*e[i-1] + e[i]

return(y)

}

plot(sim.data.armal1(100, 0.6, 0.6, 1), main=”timeseries”, xlab=”time”)

tsdisplay(sim.data.armal1(100, 0.6, 0.6, 1))

“`

#f

We make the function for the simulation of the AR(2) process

“`{r }

sim.data.ar2 <- function(n.obs, phi1, phi2, seed.nr){

set.seed(seed.nr)

y <- ts(numeric(n.obs))

e <- rnorm(n.obs)

for(i in 3: n.obs)

y[i] <- phi1*y[i-1] + phi2*y[i-2] + e[i]

return(y)

}

“`

#g

Now let’s compare processes ARIMA(1,1) and AR(2)

We see that the trends on the charts are very different. Where the last graph shows the expanding volatile dynamics.

“`{r }

plot(sim.data.armal1(100, 0.6, 0.6, 1), main=”timeseries”, xlab=”time”)

plot(sim.data.ar2(100, -0.8, 0.3, 2), main=”timeseries”, xlab=”time”)

“`

##8.6

#a-b

The function  ACF quickly decays. We see a correlation between y(t) and y(t-n). The influence of past values “y” is weakened and has a cumulative effect.

We see that the functio PACF with lag 1 goes beyond the confidence interval. And the remaining partial autocorrelations are close to zero.

The first private autocorrelation (ACF) coincides with the first autocorrelation ACF. These two bars are about the same height.

The second private autocorrelation PACF is zero. This means that yesterday’s value “y” does not carry any information about today’s values “y”.

Before us, the AR 1 process

“`{r }

df.wmurders<- wmurders

tsdisplay(df.wmurders)

“`

#c

“`{r }

model.wmurders<- arima(df.wmurders, order = c(1,0,0))

autoplot(forecast(model.wmurders))

df.wmurders.shift<- bshift(df.wmurders, k =1)

model.wmurders.shift<- arima(na.omit(df.wmurders.shift), order = c(1,0,0))

print(model.wmurders.shift)

“`

#d

We see that the remnants of both models are not particularly different. If we took a greater slip, for example 12 lags, then the remainder would be different.Also for interest I conducted a test for stationarity of the series. The test results show that we can not reject the hypothesis of stationarity of the time series.

“`{r }

checkresiduals(model.wmurders)

adf.test(df.wmurders, “stationary”)

checkresiduals(model.wmurders.shift)

“`

#e

We see that the models corrected for backshift do not differ.

“`{r }

autoplot(forecast(model.wmurders))

autoplot(forecast(model.wmurders.shift))

“`

#d

We see that the second log in ACFextends beyond the confidence interval. The histogram of the remainders is distributed not symmetrically. We can not say with certainty that this model is satisfactory.

“`{r }

checkresiduals(model.wmurders.shift)

“`

#fLet’s construct a forecast for 3 time units forward

“`{r }

autoplot(forecast(model.wmurders.shift, h = 3))

autoplot(forecast(model.wmurders, h = 3))

“`

#h

We see that the auto.arima function has a process ARIMA(1,2,1).

At the same time, the forecast graph shows a further decline. I think that this forecast is somewhat contrary to common sense, based on the previous dynamics of the time series.

“`{r }

model.wmurders.auto<- auto.arima(df.wmurders)

print(model.wmurders.auto)

print(forecast(model.wmurders.auto))

autoplot(forecast(model.wmurders.auto), h =3)

“`

##8.7

“`{r }

df.austourists<- austourists

“`

#a

We see seasonality in the time series. Also, we see an uptrend, since every new cavity of previous ones. It is also noticeable that the seasonality of the time series suffered somewhat in the period from 2001 to 2003. Most likely, the trend was influenced by some external factors.

“`{r }

autoplot(df.austourists)

tsdisplay(df.austourists)

“`

#b

We see that on the graph of the autocorrelation function each 4 lag went beyond the confidence interval. This indicates seasonality

“`{r }

Acf(df.austourists)

“`

#c

We see that on the graph of the private autocorrelation function 1, 4, 5, the lag went beyond the confidence interval. Then the function decreases sharply to zero.

We can assume that our time series is non-stationary with a season.

“`{r }

pacf(df.austourists)

“`

#d

We see that within each year in the second quarter a local minimum is formed. In addition, it is noteworthy that in 2001 and 2002 the trend was somewhat confirmed by external factors

“`{r }

ggseasonplot(df.austourists)

“`

#e

We have model ARIMA(1,0,0)(1,1,0)[4] with drift

“`{r }

model.austourists<- auto.arima(df.austourists)

print(model.austourists)

“`

#d

We see that in the first model have AIC coefficient is smaller (better) than the second model. Therefore, the first model is better.

“`{r }

model.shift<- bshift(df.austourists, k =2)

model.shift<- auto.arima(model.shift)

print(model.shift)

autoplot(forecast(model.austourists))

autoplot(forecast(model.shift))

“`

##8.8

“`{r }

df.usmelec<- usmelec

“`

#a-b

We see a waxing trend with a certain season, which has a different amplitude. At the end of the time series, we see that the seasonality amplitude is very wide and each new bottom becomes lower, the previous local bottom.

The function PACF also indicates seasonality, which occurs every 5 time units. Then the PACF function fades, which indicates that the time series is not stationary.

After the time series has been transformed, it has become more understandable. Therefore, the transformation of the time series is useful in this case.

“`{r }

autoplot(df.usmelec)

df.usmelec.ma <- SMA(df.usmelec, k = 12)

autoplot(df.usmelec.ma)

tsdisplay(df.usmelec)

tsdisplay(df.usmelec.ma)

“`

#c

It is interesting to note that for the transformed time series, we can not reject the hypothesis of stationarity.

“`{r }

Box.test(diff(df.usmelec),lag = 5, type = “Ljung-Box”)

adf.test(na.omit(df.usmelec.ma), “stationary”)

“`

#d

For the original time series, we have a process ARIMA(1,0,2)(0,1,1)[12] with drift.

For the transformed time series, we have a process ARIMA(4,0,0)(2,1,0)[12] with drift. It is important to note that for the transformed series we have the coefficient AIC 1252.7. At the same time for the original time series, we have a coefficient AIC 3050.49.

Therefore, for the transformed time series, we have the best ratio.

“`{r }

model.df.usmelec<- auto.arima(df.usmelec)

print(model.df.usmelec)

model.df.usmelec.ma <- auto.arima(df.usmelec.ma)

print(model.df.usmelec.ma)

“`

#e

In both models, the remainders are normally distributed, but the PACF goes beyond the confidence interval with some seasonality.

“`{r }

checkresiduals(model.df.usmelec)

checkresiduals(model.df.usmelec.ma)

“`

#f

We see that the model that was built on the transformed data better predicts and has smaller errors.

“`{r }

predict.model.df.usmelec<- forecast(model.df.usmelec)

autoplot(predict.model.df.usmelec)

round(accuracy(predict.model.df.usmelec))

predict.model.df.usmelec.ma <- forecast(model.df.usmelec.ma)

autoplot(predict.model.df.usmelec.ma)

round(accuracy(predict.model.df.usmelec.ma))

“`