Introduction
Moving average
Seasonal Index Method
Forecasting through simple exponential smoothing
The Holt Winters exponential smoothing forecasting model

Introduction

The following models are taken from the book “statistics for business economics / 8th edition” chapter 16 : forecasting with time series models. I am basically going through the models in the chapter one by one and implementing them in R. Some of the text / comments are coming directly from the book.

Component from a time serie to take into consideration:

Trend component
Seasonality component
Cyclical component
Irregular component

Moving average

Simple moving average

We will start with a simple moving average centered around (2m +1) on the air passenger datasets available in the datasets package. It represents the number of monthly airline passenger between 1949 and 1960.

ds <- as.integer(datasets::AirPassengers)
m <- 4

moving_avg <- function(ds){
  time_serie <- numeric(length = length(ds)) 
  
  for (j in (m+1) : (length(ds) - m)){
    X_t <- sum(ds[(j-m) : (j+m)])
    time_serie[j] <- round(1/(2*m +1) * X_t, 2)
  }
  
  return(time_serie)
  
}

time_serie <- moving_avg(ds)


ds_moving_average <- data.frame(X= 1 : length(time_serie), datasource= ds,
                                moving_avg = time_serie) 
ds_moving_average <- gather(ds_moving_average, category, value, -X)

ggplot(ds_moving_average, aes(x= X, y= value, colour= category, group= category))+
  geom_line(size= 1)+
  theme_bw()+
  theme(legend.position = "top")+
  xlab("Ratio")+
  ggtitle("Monthly Airline passagers with moving avg m= 4")

As we can see on the plot, the moving average serie is indeed, smoother than the original serie.

Weight moving average

There are several kind of moving average that can be used. For example we can usea weighted average, in which most weight is given to the central observation, with weights for other values decreasing as their distance from the central observation.

Example of a possibility for weight average with m= 2 + center and overall weight of 10

moving_avg_weight <- function(ds){
  
  time_serie <- numeric(length = length(ds))
  
  m <- 2
  
  for (j in (m+1) : (length(ds) - m)){
    X_t <- ds[(j-m) : (j+m)]
    X_t <- X_t * c(1,2,4,2,1) # multiplying by the chosen weight
    X_t <- sum(X_t)
    time_serie[j] <- round((1/10) * X_t, 2) # the overall weight is 10
  }
  return(time_serie)
}

time_serie <- moving_avg_weight(ds)

ds_moving_average_weight <- data.frame(X= 1 : length(time_serie), datasource= ds,
                                       moving_avg= time_serie)

ds_moving_average_weight <- gather(ds_moving_average_weight, category, value, -X)


ggplot(ds_moving_average_weight, aes(x= X, y= value, colour= category, group= category))+
  geom_line(size= 1)+
  theme_bw()+
  theme(legend.position = "top")+
  xlab("Ratio")+
  ggtitle("Moving avg m=2 (2*2+1) and overall weight= 10")

As there is a weight of 4 on the center, the moving average tends to follow move closely the monthly peaks

Extraction of the seasonal component through moving averages

Idea: consider a quarterly time with a seasonal component. Our strategy to remove seasonality will be to produce four periods moving averages so that the various seasonal values are brought together in a single seasonal moving average.

In the Air passenger dataset, the peaks are always during summer. So in order to remove the seasonality, we would need to use 12 rolling months.

The problem is: the location in time of the members of the series of moving averages does not correspond precisely with that of the members of the original series.

For example, the avg of x_1 x_2 x_3 and x_4 is x_2.5. So the avg of x_1 to x_12 is x_6.5

This problem can be overcome by centering our series of 11 points moving averages. This can be done by calculating the averages of adjacent pairs (x_3 = (x_2.5 + x_3.5)/2)

S <- 12 # 12 for monthly data, 4 for quarterly

moving_avg_season <- function(ds){
  
  X_t <- numeric(length = length(ds))
  time_serie <- numeric(length = length(ds))
  
  # 1. form the S point moving averages  
  for (t in (S/2) : (length(ds) - (S/2))){
    X_t[t] <- sum(ds[(t - (S/2)+1) : (t + (S/2))])/S
    # so X_t is in fact X_t + 0.5
  }
  
  # 2. form the centered S point moving average
  for (t in (S/2)+1 : (length(ds) - (S/2))){
    time_serie[t] <- (X_t[t-1] + X_t[t])/2
  }
  
  return(time_serie)
}


time_serie <- moving_avg_season(ds)

ds_moving_average_season <- data.frame(X= 1 : length(time_serie), datasource= ds, moving_avg = time_serie)

ds_moving_average_season <- gather(ds_moving_average_season, category, value, -X)


ggplot(ds_moving_average_season, aes(x= X, y= value, colour= category, group= category))+
  geom_line(size= 1)+
  theme_bw()+
  theme(legend.position = "top")+
  xlab("Ratio")+
  ggtitle("Moving avg S= 12 with seasonal adjustment")

Seasonal Index Method

Next we will discuss a seasonal adjustment approach that is based on the implicit assumption of a stable seasonal pattern over time.

We assume that for any month or quarter in each year, the effect of seasonality is to increase or decrease the series by the same percentage.

We will use our previous centered 12 points average

ds_moving_average_season <- data.frame(X= 1: length(time_serie),
                                       datasource= ds,
                                       moving_avg= time_serie)


ds_moving_average_season <- mutate(ds_moving_average_season, Percent= datasource / moving_avg * 100)

# to assess the effect of seasonality of the first month, we find the median of the 12 percentages for that month


ds_moving_average_season$Month <- rep(c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", 
                                        "Aug", "Sep", "Oct", "Nov", "Dec"), 12)


ds_moving_average_season <- ds_moving_average_season[ds_moving_average_season$Percent != Inf, ] #removing infini values

seasonal_index <- aggregate(data= ds_moving_average_season, Percent ~ Month, FUN = median)

colnames(seasonal_index) <- c("Month", "Index")

print(seasonal_index)

##    Month     Index
## 1    Apr  97.28050
## 2    Aug 120.71006
## 3    Dec  89.92974
## 4    Feb  87.37500
## 5    Jan  90.81081
## 6    Jul 125.44865
## 7    Jun 111.31915
## 8    Mar  99.54561
## 9    May  97.99692
## 10   Nov  80.19306
## 11   Oct  92.20415
## 12   Sep 105.77820

sum(seasonal_index$Index)

## [1] 1198.592

# to obtain seasonal indices, we also adjust the indices so that their average is 100

seasonal_index$Index <- seasonal_index$Index * (1200 / 1198.592)
mean(seasonal_index$Index)

## [1] 99.99999

# We interpret April figure (97) as estimating that the effect of seasonality is to April to 97% of what they would have been in the absence of seasonal factors

ds_moving_average_season <- left_join(x= ds_moving_average_season, 
                                      y= seasonal_index, by.x= Month)

## Joining, by = "Month"

# finally we obtain our seasonality adjusted value = original value * (100/Index)

ds_moving_average_season <- mutate(ds_moving_average_season, adjusted_value = datasource * (100 / Index))


head(ds_moving_average_season)

##    X datasource moving_avg   Percent Month     Index adjusted_value
## 1  7        148   126.7917 116.72691   Jul 125.59602       117.8381
## 2  8        148   127.2500 116.30648   Aug 120.85186       122.4640
## 3  9        136   127.9583 106.28460   Sep 105.90246       128.4201
## 4 10        119   128.5833  92.54699   Oct  92.31246       128.9100
## 5 11        104   129.0000  80.62016   Nov  80.28726       129.5349
## 6 12        118   129.7500  90.94412   Dec  90.03538       131.0596

plot(ds_moving_average_season$X, ds_moving_average_season$adjusted_value,
     xlab = "Timeline", 
     ylab = "Adjusted Serie",
     main= "Seasonality adjusted air passenger numbers")

Forecasting through simple exponential smoothing

Exponential smoothing is appropriate when the series is non seasonal and has no consistent upward or downward trend.

Simple exponential smoothing allows a compromise between these extremes, providing a forecast based on weighted average of current and past values. In forming this average, most weight is given to the most recent observation, rather less to the immediate preceding value, less to the the one before that, and so on.

Note that a small alpha value gives greater weight to the forecast x_t-1, which is based on the past history of the series

# we use a new dataset
plot(datasets::BJsales)

ds <- data.frame(X = 1: length(datasets::BJsales), datasource= datasets::BJsales, alpha06= 0, alpha01 = 0)  

# example with alpha = 0.6

ds$alpha06[1] <- ds$datasource[1] # first forecast value is equal to the first value

alpha <- 0.6

for (i in 2: length(datasets::BJsales)){
  ds$alpha06[i] <- (1-alpha) * ds$alpha06[i - 1] + (alpha * ds$datasource[i]) 
}

# example with alpha = 0.1

ds$alpha01[1] <- ds$datasource[1]
alpha <- 0.1

for (i in 2: length(datasets::BJsales)){
  ds$alpha01[i] <- (1-alpha) * ds$alpha01[i - 1] + (alpha * ds$datasource[i]) 
}

print(head(ds))

##   X datasource  alpha06  alpha01
## 1 1      200.1 200.1000 200.1000
## 2 2      199.5 199.7400 200.0400
## 3 3      199.4 199.5360 199.9760
## 4 4      198.9 199.1544 199.8684
## 5 5      199.0 199.0618 199.7816
## 6 6      200.2 199.7447 199.8234

#plot

time_serie <- gather(ds, category, value, -X)

## Warning: attributes are not identical across measure variables;
## they will be dropped

ggplot(time_serie, aes(x= X, y= value, colour= category, group= category)) + 
  geom_line(size= 1)+
  theme_bw()+
  theme(legend.position = "top")+
  xlab("Time")+
  ggtitle("Smoothing plot for BJsales for alpha = 0.1 and 0.6")

The Holt Winters exponential smoothing forecasting model

The Holt Winters exponential smoothing procedure allows for trend, and possibly also seasonality, in a time serie.

Holt Winters exponential smoothing using R function

Note: code found on https://www.r-bloggers.com/holt-winters-forecast-using-ggplot2/

demand <- ts(datasets::BJsales, start = c(2000, 1), frequency = 12)

hw <- HoltWinters(demand)

plot(hw)

# next we calculate the forecast for 12 months with a confidence interval of 0.95 and plot the forecast together with actual and fitted values.

forecast <- predict(hw, n.ahead = 12, prediction.interval = T, level = 0.95)

plot(hw, forecast)

Holt Winters exponential smoothing

The principle is similar to the simple exponential smoothing. We will use 2 smoothing constants, Alpha and Beta (between 0 and 1) as well as 2 equations:

Forecast of x_t
Trend estimate represented by T_t

Those equations are used to update previous estimate using new observations

ds <- data.frame(X= 1: length(datasets::BJsales), datasource= datasets::BJsales, forecast= 0, trend= 0)

# initial values
ds$forecast[2] <- ds$datasource[2]
ds$trend[2] <- ds$datasource[2] - ds$datasource[1]

# Obtain estimate of level forecast X_t and trend T
alpha <- 0.3
beta <- 0.3

for (i in 3 : length(datasets::BJsales)){
  
  ds$forecast[i] <- (1 - alpha) * (ds$forecast[i - 1]) + ds$trend[i -1] + 
    alpha * ds$datasource[i]
  
  ds$trend[i] <- (1 - beta) * ds$trend[i - 1] + beta * (ds$forecast[i] - ds$forecast[i - 1])
  
} 

# The forecast is obtained by forecast x_n+h = forecast x_n + h * T_n

# We are standing at a time n and h is the number of period in the future

forecast_serie <- data.frame(X = 151:161, forecast= 0)


for (i in 1:11){
  
  forecast_serie$forecast[i] <- tail(ds$forecast, 1) + i * tail(ds$trend, 1)
  
}

ds <- bind_rows(ds, forecast_serie)

## Warning in bind_rows_(x, .id): Vectorizing 'ts' elements may not preserve
## their attributes

# plot

ds <- ds[, 1:3] # drop trend

time_serie <- gather(ds, category, value, -X)

ggplot(time_serie, aes(x= X, y= value, colour= category, group= category))+
  geom_line(size= 1)+
  theme_bw()+
  theme(legend.position = "top")+
  xlab("Time")+
  ggtitle("Holt Winters for alpha = 0.3 and beta = 0.3")+
  scale_y_continuous(limits = c(180, max(time_serie$value) + 10))

## Warning: Removed 12 rows containing missing values (geom_path).

Holt Winters exponential smoothing with seasonal factor

3 constants are used between 0 and 1: alpha, beta and gamma. The closer a parameter is to 0, the lower the weight of the present observation and the higher the weight of previous estimates in determining the updated statistic.

Gamma is the multiplicative seasonal factor. S is used to denote period. So 4 for quarterly data and 12 for monthly data.

ds <- data.frame(X= 1 : length(datasets::BJsales), datasource= datasets::BJsales, forecast= 0, trend= 0, season= 0)

alpha <- 0.7
beta <- 0.6
gamma <- 0.3
S <- 12


# intitial values
# There are different ways to initiate values. We will use the one found on <https://www.itl.nist.gov/div898/handbook/pmc/section4/pmc435.htm>

# intial value for x forecast
# because we have to do a full year for the season, we will start at the 12th value
# as for previous Holt Winter method we will set it up with the X value

ds$forecast[12] <- ds$datasource[12]

# initial value for Trend factor
ds$trend[12] <- (1/S) * sum((ds$datasource[(S+1) : (S+S)] - ds$datasource[1:12])/S)

# initial value for seasonal factor
# step 1: compute average of each year 12 years in our case as the 13th year is not complete

seasonal_fct <- numeric(length = 12)

for (i in 1:12){
  seasonal_fct[i] <- mean(ds$datasource[(i * 12 -11): (i*12)])
}

# step 2: divide observation by the appropriate yearly mean
seasonal_fct <- rep(seasonal_fct, each= 12, times= 1)
seasonal_fct <- ds$datasource[1:144] / seasonal_fct

# step 3
seasonal_fct <- matrix(data = seasonal_fct, ncol = 12, byrow = FALSE)
seasonal_fct <- rowSums(seasonal_fct) / 12

ds$season[1:12] <- seasonal_fct

# calculate forecast, trend and seasonal

for (i in 13:150){
  
  ds$forecast[i] <- (1 - alpha) * (ds$forecast[i - 1] + ds$trend[i - 1]) + alpha * (ds$datasource[i] / ds$season[i - S])
  
  ds$trend[i] <- (1 - beta) * ds$trend[i - 1] + beta * (ds$forecast[i] / ds$forecast[i-1])
  
  ds$season[i] <- (1 - gamma) * ds$season[i - S] + gamma * (ds$datasource[i] / ds$forecast[i])
}

# the forecast is obtained by forecast x_n+h = (forecast x_n + hT_n) * F_n+h-S
# F is the one generated for the most recent time period 
# we forecast future values at h time periods ahead of the last observation x_n

forecast_serie <- data.frame(X= 151:161, forecast= 0)

for (i in 1*11){
  
  forecast_serie$forecast[i] <- (tail(ds$forecast, 1) + i * tail(ds$trend, 1)) * ds$season[150 + i - S]
  
}

ds <- bind_rows(ds, forecast_serie)

## Warning in bind_rows_(x, .id): Vectorizing 'ts' elements may not preserve
## their attributes

# plot
ds <- ds[-c(1:12), -c(4,5)] # drop tgrend and seasonal

time_serie <- gather(ds, category, value, -X)

ggplot(time_serie, aes(x= X, y= value, colour= category, group= category))+
  geom_line(size= 1)+
  theme_bw()+
  theme(legend.position = "top")+
  xlab("Time")+
  ggtitle("Holt Winters for Alpha= 0.7, Beta= 0.6, Gamma= 0.3")+
  scale_y_continuous(limits = c(180, max(time_serie$value) + 10))

## Warning: Removed 11 rows containing missing values (geom_path).

The alpha, beta and gamma values strongly impact the forecast. To find the best combination, the sum of squared errors casn be used

Forecasting with Time serie models

ndee

May 28, 2018