---
title: "Practice / Forecasting"
author: "Martin Schedlbauer"
output:
  html_document:
    df_print: paged
---

## Task Set I

### Question I.1

The built-in dataset [USArrests](https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/USArrests.html) contains statistics about violent crime rates in the US States.

Let's start by inspecting the data set's structure.

```{r}
head(USArrests)
```

We want to determine which states are outliers in terms of murders. Outliers, for the sake of this question, are defined as values that are more than 1.5 standard deviations from the mean.

```{r}
df <- USArrests

# calculate mean and standard deviation
mean.murders <- mean(df$Murder)
sd.murders <- sd(df$Murder)

# find rows that contains outliers
outliers <- which(abs((df$Murder - mean.murders) / sd.murders) > 1.5)

# the states are not actually columns, but are row names
# row names are a unique concept for R data frames
rownames(df[outliers,])
```

## Question I.2

For the same dataset, is there a correlation between urban population and murder, *i.e.*, as one goes up, does the other statistic as well? Comment on the strength of the correlation. Which correlation algorithm is appropriate? Pearson? Spearman, Kendall? How would you decide between them? What if you choose an incorrect algorithm; what would the effect be?

The most common way to calculate correlation is the Pearson Coefficient of Correlation. However, this metric presumes that the data is (reasonably) normally distributed. We can either run a normality test such as the Shapiro-Wilk test or plot the data in a histogram or in a QQ plot and visually inspect the data for normality.

Let's start with a QQ plot. A QQ plot (or quantile-quantile plot) draws the correlation between a given sample and the normal distribution. A 45-degree reference line is also plotted. QQ plots are used to visually check the normality of the data. The data should closely follow the 45-degree line if it is normally distributed.

```{r}
qqnorm(df$Murder, pch = 1, frame = FALSE)
qqline(df$Murder, col = "steelblue", lwd = 2)

```

The plot shows some outliers, so let's run a Shapiro-Wilk test as well. For it to be normally distributed, its *W* should be close to 1 and the p-value \> 0.05.

```{r}
shapiro.test(df$Murder)
```

From the QQ Plot and the Shapiro-Wilk test, we can conclude that the data is reasonably normally distributed and thus the Pearson correlation metric is appropriate to calculate. Note that one could always calculate the metric, but if the data is not normally distributed then its result could lead us to falsely conclude that there is a correlation.

```{r}
cor(df$Murder, df$UrbanPop, method = "pearson")
```

For there to be a correlation, the coefficient must be close to 1. In our case, it is close to 0, so there is no correlation between murders and urban population.

## Task Set II

### Question II.1

Based on the [data on the growth of mobile phone use in Brazil](https://docs.google.com/spreadsheets/d/1tOnM9XceK4Ak8tzWQ2vDelWlJexzJiS3LbT6MN6_rW0/edit?usp=sharing) (you'll need to copy the data and create a CSV that you can load into R or use the `gsheet2tbl()` function from the **gsheet** package), forecast phone use for the next time period using a 2-year weighted moving average (with weights of 5 for the most recent year, and 2 for other), exponential smoothing (alpha of 0.4), and linear regression trendline.

```{r warning=F, message=F}
# be sure to install the package if not installed
library(gsheet)

url <- "https://docs.google.com/spreadsheets/d/1tOnM9XceK4Ak8tzWQ2vDelWlJexzJiS3LbT6MN6_rW0/edit?usp=sharing"
df <- read.csv(text=gsheet2text(url, format='csv'), 
               stringsAsFactors=FALSE,
               encoding = "UTF-8")
```

Alternatively, you could have opened the data in a Google Sheet by following the URL and then download the data as a CSV; and then load the local CSV file (after moving it to the project folder, of course).

Let's inspect the data.

```{r}
head(df,4)
tail(df,3)
```

Notice there is a 12th row -- that is the time period for which we need a forecast, so let's delete that from the data. In R, there is no actual way to "delete" a row, so we copy the data to a new dataframe and omit the last row from the copy process. Notice how we assign the 11 rows from the old dataframe to a variable that's the same as the old one -- we are reusing and redefining -- this is common in programming, but perhaps not always the best choice as it overwrites the old dataframe. For large dataframes this has the added benefit of releasing the memory for the old dataframe -- and be cognizant of memory use for large data sets is critical.

```{r}
df <- df[1:(nrow(df)-1),]

tail(df,3)
```

### Weighted Moving Average

Now that we have the data we can start building our forecasts. Start with the simplest of the models (weighted average) and build the code first. Then we'll make the code more generic and more reusable by moving the forecast calculation into a function. Of course we will have to test the code and the function to make sure they work as intended and make a forecast.

A weighted moving average takes the average of some prior observations. We will forecast phone use for the next time period using a 2-year weighted moving average (with weights of 5 for the most recent year, and 2 for other). Note that for a weighted average we need to multiply the values by their weights and divide by the sum of the weights.

```{r}
n <- nrow(df)
w1 <- 5
w2 <- 2
f <- ((w1 * df$Subscribers[n]) + (w2 * df$Subscribers[n-1])) / (w1 + w2)

print(f)
```

Using a 2-year weighted average with weights of `r w1` and `r w2`, we get a forecast of `r format(f, big.mark = ",", scientific = F)` for the next year.

Let's turn this into a function -- let's call the function `wmaForecast()`. We will start by moving the code into a function and taking the time series data (as a vector) and the weights as arguments (parameters) for the function so that we have configurability.

```{r}
wmaForecast <- function (data, w1, w2) {
  f <- ((w1 * data[n-1]) + (w2 * data[n-2])) / (w1 + w2)
  return (f)
} 
```

Next, we'll test the function by calling it with the same weights and we, of course, get the same answer.

```{r}
next.year.wma <- wmaForecast(df$Subscribers, 5, 2)

print(next.year.wma)
```

Now let's make the function even more generic by allowing us to specify the number of time periods and get the weights in a vector.

```{r}
wmaForecast <- function (data, n, w) {
  
  # get the indexes for the periods
  periods <- length(data):(length(data) - n + 1)
  
  # multiply the periods by the weights and sum
  p <- sum(data[periods] * w)
  
  # divide the sum by the sum of the weights
  f <- p / sum(w)
  
  # return forecast
  return (f)
} 
```

Let's first test this with the same values as before to ensure we get the same answer.

```{r}
next.year.wma <- wmaForecast(df$Subscribers, 2, c(5,2))
print(next.year.wma)
```

Now let's test this with a three-year moving average and weights of (5,3,2). You should verify that the answer is correct by manually calculating the forecast (I did that and it corresponds to the value returned by the function).

```{r}
next.year.3y <- wmaForecast(df$Subscribers, 3, c(5,3,2))
print(next.year.3y)
```

### Exponential Smoothing

This method of making a forecast required bootstrapping and it is recursive which means that the forecast depends on the prior forecast which, in turn, depends on the prior until the beginning. So, we need to start with the first time period for which we can calculate a "forecast" and bootstrap our way forward. We will do this "brute-force" but a more elegant programming solution would be to build a recursive function.

Let's start again by writing the code to calculate a forecast and then turn that code into a function we can call.

```{r}
# exponential smoothing starts at the first period for which we have a forecast
# we need to calculate the error of the forecast against the actual observed value

# the "forecast" for the first time period is the observed value (by convention)
df$f[1] <- df$Subscribers[1]
# the error is (by definition) then 0
df$error[1] <- 0

alpha <- 0.4

# calculate the forecast for all period > 1
for (i in 2:(nrow(df))) {
  # calculate the forecast as forecast(t-1) + alpha*error(i-1)
  df$f[i] <- df$f[i-1] + (alpha * df$error[i-1])
  
  # calculate the error for the forecast (Y - f)
  df$error[i] <- df$Subscribers[i] - df$f[i]
}

next.year.es <- df$f[i] + (alpha * df$error[i])
print(next.year.es)
```

We got a forecast; it's a bit lower than the other two, but it is "in line". So, now we can convert out code into a more generic function. Once again, we expect a vector of time-series data as as argument, plus a value for "alpha".

```{r}
esForecast <- function (data, alpha) {
  # we need the data in a dataframe so we can have columns for forecast (f) and error (e)
  d.es <- data.frame(t = 1:length(data),
                     x = data,
                     f = 0,
                     e = 0)
  
  # the "forecast" for the first time period is the observed value (by convention)
  d.es$f[1] <- d.es$x[1]
  # the error is (by definition) then 0
  d.es$e[1] <- 0
  
  # calculate the forecast for all period > 1
  for (i in 2:(length(data))) {
    # calculate the forecast as forecast(t-1) + alpha*error(i-1)
    d.es$f[i] <- d.es$f[i-1] + (alpha * d.es$e[i-1])
    
    # calculate the error for the forecast (Y - f)
    d.es$e[i] <- d.es$x[i] - d.es$f[i]
  }
  
  return (d.es$f[i] + (alpha * d.es$e[i]))
}
```

To test, we will call the function with our data:

```{r}
next.year.es <- esForecast(df$Subscribers, 0.4)
print(next.year.es)
```

As before, same value... so we have at least some confidence that our function "works".

### Trend Line

The last forecasting method we will apply is a linear regression trendline. We will use the `lm()` function of R, although other functions could also be used.

```{r}
forecast.model <- lm(Subscribers ~ Year, 
                     data = df)
```

Now that we have a "linear model", we can access the intercept and slope and build a linear equation in the form ($y = mx + b$) and calculate the forecast by using *x = 12*, the "next year".

```{r}
f.lm <- forecast.model$coefficients[[2]] * 12 + forecast.model$coefficients[[1]]
print(f.lm)
```

Note that the model has an element called "coefficients" which is a *list* of coefficients (intercept and slope), so we need to access them with the list access operator [[]].

The forecast using linear regression trend line is close to what we got for the other two. Of course, once again, we need to create a function. In the function, we will set "time" to a sequence starting at 1, so the data we get is just a vector of prior time period values.

```{r}
tlForecast <- function (data, t) {
  periods = 1:length(data)
  
  # lm() requires them to be in a dataframe
  d.lm <- data.frame(x = periods,
                     y = data)
  
  m <- lm(y ~ x, data = d.lm)
  
  f.lm <- m$coefficients[[2]] * (t) + m$coefficients[[1]]
  
  return (f.lm)
}
```

And, as before, let's compare to what we got before. As you can see it matches what we got before (as it should).

```{r}
next.year.tl <- tlForecast(df$Subscribers, 12)
print(next.year.tl)
```

Now, we have three ways to calculate a forecast using three different methods. Later we will define a function that combines the forecast of the three methods through a weighted average. That way we get a more accurate forecast.

## Question II.2

Calculate the squared error for each model, *i.e.*, use the model to calculate a forecast for each provided time period in the data set and then the square the error.

### Model Validation

A common approach to "testing" a model is to apply the model to the data that we have and see how well the model's prediction matches the actual data. This means we apply the model to every value we have in the time-series and for which a prediction is possible. For example, since the WMA requires 2 years (or whichever the time periods are), so we can only calculate a "forecast" for years 3 and later. Once we have a prediction, we calculate the squared error (the square of the difference between the actual and the predicted value), *MSE*. Other metrics of model fit is the mean absolute error (or deviation), or *MAD*.

```{r}
# calculate the MSE for each model

# initialize new columns for MSE (this is important; R thing)
df$mse.TL <- 0
df$mse.ES <- 0
df$mse.MA <- 0

## calculate MSE for E/S: need to start at 3 as 1 is the "bootstrap"
for (i in 3:nrow(df)) {
  # forecast ith value with E/S
  f.ES <- esForecast(df$Subscribers[1:(i-1)], alpha = 0.4)
  
  # calculate MSE
  df$mse.ES[i] <- (df$Subscribers[i] - f.ES)^2
}

mse.ES <- mean(df$mse.ES)
print(paste0("MSE for E/S: ", mse.ES))

## calculate MSE for T/L
for (i in 1:nrow(df)) {
  # forecast ith value with T/L
  f.TL <- tlForecast(df$Subscribers, i)
  
  # calculate MSE
  df$mse.TL[i] <- (df$Subscribers[i] - f.TL)^2
}

mse.TL <- mean(df$mse.TL)
print(paste0("MSE for TL: ", mse.TL))

## calculate MSE for WMA: need to start at 3
for (i in 3:nrow(df)) {
  # forecast ith value with WMA
  f.WMA <- wmaForecast(df$Subscribers[1:(i-1)], 2, c(5,2))
  
  # calculate MSE
  df$mse.MA[i] <- (df$Subscribers[i] - f.WMA)^2
}

mse.MA <- mean(df$mse.MA)
print(paste0("MSE for MA: ", mse.MA))
```

As an aside: how do we find the "optimal" *alpha* or the optimal weights and number of time periods to average in WMA? Answer: simple, we try lots of different alpha's; for each we calculate the MSE and/or MAD and we choose the alpha that results in the smallest MSE/MAD. Same for the weights: we use our computers to find the "optimal" solution: machine learning. Of course, this "training" and parameter calibration takes time, especially for large data sets.

From the results, it appears that the Trend Line model has the smallest MSE, so it most closely matches its historical data, *i.e.*, it has the best "backfit". So, the predictions from this model are likely to be most accurate.

### Ensemble Forecast

We will write a function called `ensembleForecast()` that calculates a weighted average forecast by averaging out the three forecasts calculated with the following weights: 4 for trend line, 2 for exponential smoothing, 1 for weighted moving average. The function takes a value for *alpha* and the weights and number of time periods for the weighted average forecast (*weights.WMA* and *n*) as well as the weights for the ensemble (*w*).

```{r}
ensembleForecast <- function (data, t, alpha, n, weights.WMA, w) {
  f.wma <- wmaForecast(data, n, weights.WMA)
  f.es <- esForecast(data, alpha)
  f.tl <- tlForecast(data, t)
  
  f <- (w[1]*f.tl + w[2]*f.es + w[3]*f.wma) / sum(w)
}
```

Now that we have the ensemble function, we can call it and get a forecast that combines three separate forecasts.

```{r}
next.year.ensemble <- ensembleForecast(df$Subscribers, 12, 0.4, 2, c(5,2), c(4,2,1))
print(next.year.ensemble)
```
