R Lab - Time Series
Tanvir
2023-12-11
rm(list = ls())R is very powerful when it comes to Time Series Analysis, however the methods that we are looking at this chapter are very basic, so probably you won’t appreciate the use of R for these methods. Sorry guys…. ECO204 is just the starting, so we can’t cover the advanced methods here and also I think some of the concepts of time series is beyond the scope of this course. However, if you are interested in Time Series Analysis, then you can take some courses in Econometrics and hopefully you will learn more. Nevertheless, here are some resources,
Time Series Objects and Plots
In R for time series analysis there is a special object. We can
convert any vector to time series objects by calling the function
ts(). The ts() function takes two arguments,
the first one is the data and the second one is the frequency of the
data. If we have consecutive time periods for a short amount of time,
probably you don’t need to do anything with the frequency, but
complicated time structure like seasonality or cycles, possibly you
should be careful.
Okay let’s see an example, where we will start with a vector and then convert it to a time series object.
raw_data <- c(17, 21, 19, 23, 18, 16, 20, 18, 22, 20, 15, 22, 31, 34, 31, 33, 28, 32, 30, 29, 34, 33)
# here frequency is 1 since we observe the data in every week, we don't need to
# set this but I did it just for clarity
gas_revised <- ts(raw_data, frequency = 1)
# "o" is used to show the circle in the plots
plot.ts(gas_revised , xlim = c(0, 22), ylim = c(0, 40), type = "o",
xlab = "Week", ylab = "Sales")
One very important thing in the time series is the lagged values. We
can get the lagged values of a time series object with the function
lag(). Careful this function will only work for time series
objects not for any arbitrary vector. Let’s see how we can use this
function.
# This will creae a new object called gas_lag which is the lagged value of the gas_revised object
lag(gas_revised, k = -1)## Time Series:
## Start = 2
## End = 23
## Frequency = 1
## [1] 17 21 19 23 18 16 20 18 22 20 15 22 31 34 31 33 28 32 30 29 34 33Note, we used k = -1. The way you should understand
this, we are in the period “0” currently. Since we want to get the data
for the past period, so we used k = - 1. If we use
k = 1 then we will see value in the future, i.e., one
period ahead. We can use the function cbind() to combine
the original time series object and the lagged time series data.
cbind(gas_revised, lag(gas_revised, k = -1))## Time Series:
## Start = 1
## End = 23
## Frequency = 1
## gas_revised lag(gas_revised, k = -1)
## 1 17 NA
## 2 21 17
## 3 19 21
## 4 23 19
## 5 18 23
## 6 16 18
## 7 20 16
## 8 18 20
## 9 22 18
## 10 20 22
## 11 15 20
## 12 22 15
## 13 31 22
## 14 34 31
## 15 31 34
## 16 33 31
## 17 28 33
## 18 32 28
## 19 30 32
## 20 29 30
## 21 34 29
## 22 33 34
## 23 NA 33We can also different lagged values, for different k, you can try
this. In R we can also look at directly the difference between the time
series object and the lagged time series object. We can do that using
the diff() function. Since I don’t need this I won’t use
this function now, maybe later!
Horizontal Pattern
When the plot shows a horizontal pattern, then we can say that probably the time series object is stationary. For Horizontal pattern we learned two techniques to forecast,
- Naive Forecast
- Moving Average Forecast (Equal or Unequal Weights)
Let’s see them now
1. Naive Forecast
In this method we simply take the last observation (lagged observation) as the forecast for the current period. So the equation is
\[ F_t = y_{t - 1} \]
This means that the forecast for the current period is the last observation. Let’s see how we can do this in R. We will use the gasoline data from the Anderson book. First we will load the data, create a time series object and plot it
library(readxl)
raw_data <- read_excel("/home/tanvir/Documents/ownCloud/Git_Repos/EWU_repos/3_Fall_2023/eco_204/ewu-eco204.github.io/pdf_files/slides/06_time_series/Codes/Gasoline.xlsx")
# take only the sales column and convert it to a time series object
gas_data <- ts(raw_data$Sales, frequency = 1)
# see the data
gas_data## Time Series:
## Start = 1
## End = 12
## Frequency = 1
## [1] 17 21 19 23 18 16 20 18 22 20 15 22Let’s plot the time series object
plot.ts(gas_data, type = "o", xlab = "Week", ylab = "Sales")
Now we will do the forecast using the naive method and then see
MFE, MAE or other summary measures of the
forecast. We will use the accuracy() function from the
package forecast to do that. If you don’t have the package
installed, then you can install it using the command
install.packages("forecast")
# get the naive forecast
gas_forecast_naive <- lag(gas_data, k = -1)
# load the library forecast
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo# get the summary measures for accuracy of the forecast
accuracy(gas_forecast_naive, gas_data)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.4545455 4.033947 3.727273 0.1082169 19.24431 -0.4553872 1Note that the accuracy() function takes two arguments,
the first one is the forecast and the second one is the actual data. The
output of the accuracy() function is a table, which shows
the summary statistics of the forecast.
The first column
MEis Mean Forecast Error (MFE), which is the average of the forecast errors.The second column
RMSEis Root Mean Square Error (RMSE), which is the square root of the average of the squared forecast errors. If you square the RMSE, then you will get the MSE (Mean Squared Error).The third column
MAEis Mean Absolute Error (MAE), which is the average of the absolute forecast errors.The fourth column
MPEis Mean Percentage Error (MPE), which is the average of the percentage forecast errors.The fifth column
MAPEis Mean Absolute Percentage Error (MAPE), which is the average of the absolute percentage forecast errors.
You can ignore the last two columns. Finally we can plot the forecast
and the actual data together to see how the forecast looks like. In this
case we will use the plot.ts() function. To do this first
we need to make a matrix of the forecast and the actual data using the
cbind() function. Then we will use the
plot.ts() function to plot the data. Let’s see how we can
do that.
# make a matrix of the forecast and the actual data
gas_forecast_naive_m <- cbind( gas_data, gas_forecast_naive)
gas_forecast_naive_m## Time Series:
## Start = 1
## End = 13
## Frequency = 1
## gas_data gas_forecast_naive
## 1 17 NA
## 2 21 17
## 3 19 21
## 4 23 19
## 5 18 23
## 6 16 18
## 7 20 16
## 8 18 20
## 9 22 18
## 10 20 22
## 11 15 20
## 12 22 15
## 13 NA 22# remove the NAs from the matrix
gas_forecast_naive_m <- na.omit(gas_forecast_naive_m)
gas_forecast_naive_m ## Time Series:
## Start = 2
## End = 12
## Frequency = 1
## gas_data gas_forecast_naive
## 2 21 17
## 3 19 21
## 4 23 19
## 5 18 23
## 6 16 18
## 7 20 16
## 8 18 20
## 9 22 18
## 10 20 22
## 11 15 20
## 12 22 15# plot the data
plot.ts(gas_forecast_naive_m, plot.type = "single", xlab = "Week", ylab = "Sales", col = c("red", "blue"), main = "Naive Forecast", type = "o", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))
2. Moving Average Forecasts (Equal / Unequal Weights)
In this method we take the average of the last few observations as the forecast for the current period. The number of observations that we take the average is called the order. For example, if we take the average of the last 3 observations, then the order is 3. If we write this then,
\[ F_t = \frac{y_{t - 1} + y_{t - 2} + y_{t - 3}}{3} \]
We will use the function rollmean() from the package
xts to do that. Let’s see how we can do that.
library(xts)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericgas_ma <- rollmean(gas_data, k = 3, align = "right", fill = NA)
gas_ma## Time Series:
## Start = 1
## End = 12
## Frequency = 1
## [1] NA NA 19 21 20 19 18 18 20 20 19 19Now this is the data set of forecasts. The problem with the
rollmean() function is that it does not give us 3 NAs at
the beginning of the data set (which we should have), it gives two, but
we can fix that using the lag() function. So using
lag and then cbind() we can get the data set
of forecasts.
cbind(gas_data, lag(gas_ma, k = -1))## Time Series:
## Start = 1
## End = 13
## Frequency = 1
## gas_data lag(gas_ma, k = -1)
## 1 17 NA
## 2 21 NA
## 3 19 NA
## 4 23 19
## 5 18 21
## 6 16 20
## 7 20 19
## 8 18 18
## 9 22 18
## 10 20 20
## 11 15 20
## 12 22 19
## 13 NA 19Finally we can use the accuracy() function to calculate
the summary statistics of the forecast.
# fix the lag
gas_forecast_ma <- lag(gas_ma, k = -1)
# get the summary measures for accuracy of the forecast
accuracy(gas_forecast_ma, gas_data)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0 3.197221 2.666667 -2.310057 14.35661 -0.2065217 0.7135559Also here we can plot the forecast and the actual data together to see how the forecast looks like.
# make a matrix of the forecast and the actual data
gas_forecast_ma_m <- cbind( gas_data, gas_forecast_ma)
# remove the NAs from the matrix
gas_forecast_ma_m <- na.omit(gas_forecast_ma_m)
gas_forecast_ma_m ## Time Series:
## Start = 4
## End = 12
## Frequency = 1
## gas_data gas_forecast_ma
## 4 23 19
## 5 18 21
## 6 16 20
## 7 20 19
## 8 18 18
## 9 22 18
## 10 20 20
## 11 15 20
## 12 22 19# plot the data
plot.ts(gas_forecast_ma_m, plot.type = "single", xlab = "Week", ylab = "Sales", col = c("red", "blue"), main = "MA Forecast (Equal Weights)", type = "o", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))
Finally we will do the moving average with unequal weights. For this we can just give weights for the three observations that we are taking the average. Let’s see how we can do that.
gas_ma_unequal <- rollapply(gas_data, width = 3, FUN = function(x){weighted.mean(x, w = c(1/6, 2/6, 3/6))}, align = "right", fill = NA)
gas_ma_unequal## Time Series:
## Start = 1
## End = 12
## Frequency = 1
## [1] NA NA 19.33333 21.33333 19.83333 17.83333 18.33333 18.33333
## [9] 20.33333 20.33333 17.83333 19.33333Again we can use the accuracy() function to calculate
the summary statistics of the forecast. Let’s see how we can do
that.
accuracy(lag(gas_ma_unequal, k = -1), gas_data)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0.05555556 3.3898 2.981481 -2.129945 15.99248 -0.2994209 0.7723224Again Finally we can plot the forecast and the actual data together to see how the forecast looks like.
# make a matrix of the forecast and the actual data
gas_forecast_ma_unequal_m <- cbind( gas_data, lag(gas_ma_unequal, k = -1))
# remove the NAs from the matrix
gas_forecast_ma_unequal_m <- na.omit(gas_forecast_ma_unequal_m)
gas_forecast_ma_unequal_m## Time Series:
## Start = 4
## End = 12
## Frequency = 1
## gas_data lag(gas_ma_unequal, k = -1)
## 4 23 19.33333
## 5 18 21.33333
## 6 16 19.83333
## 7 20 17.83333
## 8 18 18.33333
## 9 22 18.33333
## 10 20 20.33333
## 11 15 20.33333
## 12 22 17.83333# plot the data
plot.ts(gas_forecast_ma_unequal_m, plot.type = "single", xlab = "Week", ylab = "Sales", col = c("red", "blue"), main = "MA Forecast (Unequal Weights)", type = "o", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))
3. Exponential Smoothing
We didn’t cover exponential smoothing!
Trend Pattern
Now we will forecasting with trend pattern. For trend pattern we can use forecast using linear or non-linear trend pattern.
1. Linear Trend Pattern
Let’s load the Bicycle data from the Anderson book. As usual, first we will load the data, create a time series object and plot it
raw_data <- read_excel("/home/tanvir/Documents/ownCloud/Git_Repos/EWU_repos/3_Fall_2023/eco_204/ewu-eco204.github.io/pdf_files/slides/06_time_series/Codes/Bicycle.xlsx")
bicycle_data <- ts(raw_data$Sales, frequency = 1)
plot.ts(bicycle_data, type = "o", xlab = "Year", ylab = "Sales")
The data shows a trend pattern. Now we will use the
tslm() function to fit a linear trend model. Let’s see how
we can do that.
bicycle_lm <- tslm(bicycle_data ~ trend)
summary(bicycle_lm)##
## Call:
## tslm(formula = bicycle_data ~ trend)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.900 -1.275 0.200 0.500 3.400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 20.4000 1.3382 15.24 3.4e-07 ***
## trend 1.1000 0.2157 5.10 0.00093 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.959 on 8 degrees of freedom
## Multiple R-squared: 0.7648, Adjusted R-squared: 0.7354
## F-statistic: 26.01 on 1 and 8 DF, p-value: 0.0009295Now we will use the predict() function to forecast the
data. Let’s see how we can do that.
bicycle_data_forecast_table <- forecast(bicycle_lm)
bicycle_data_forecast_table## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 11 32.5 29.18618 35.81382 27.02921 37.97079
## 12 33.6 30.12574 37.07426 27.86434 39.33566
## 13 34.7 31.04746 38.35254 28.67002 40.72998
## 14 35.8 31.95382 39.64618 29.45033 42.14967
## 15 36.9 32.84701 40.95299 30.20891 43.59109
## 16 38.0 33.72895 42.27105 30.94893 45.05107
## 17 39.1 34.60128 43.59872 31.67306 46.52694
## 18 40.2 35.46538 44.93462 32.38362 48.01638
## 19 41.3 36.32243 46.27757 33.08252 49.51748
## 20 42.4 37.17340 47.62660 33.77140 51.02860# we only need the point forecast now
bicycle_data_forecast <- bicycle_data_forecast_table$fittedNow we will calculate the summary statistics of the forecast. Same technique as before.
accuracy(bicycle_data_forecast, bicycle_data)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -3.552844e-16 1.752142 1.32 -0.4427101 5.068143 0.08794788 0.6752526Finally we will plot the actual data and the forecasted data together.
bicycle_data_forecast_m <- cbind(bicycle_data, bicycle_data_forecast)
bicycle_data_forecast_m <- na.omit(bicycle_data_forecast_m)
plot.ts(bicycle_data_forecast_m , type = "o", plot.type = "single", xlab = "Year", ylab = "Sales", col = c("red", "blue"), main = "Linear Trend Forecast", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))
Finally we can fit a quadratic trend model. Let’s see how we can do that. In this case we can use the Cholesterol data from the Anderson book.
raw_data <- read_excel("/home/tanvir/Documents/ownCloud/Git_Repos/EWU_repos/3_Fall_2023/eco_204/ewu-eco204.github.io/pdf_files/slides/06_time_series/Codes/Cholesterol.xlsx")
cholesterol_data <- ts(raw_data$Revenue , frequency = 1)
plot.ts(cholesterol_data, type = "o", xlab = "Year", ylab = "Cholesterol")
cholesterol_lm <- tslm(cholesterol_data ~ trend + I(trend^2))
summary(cholesterol_lm)##
## Call:
## tslm(formula = cholesterol_data ~ trend + I(trend^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6767 -2.2459 -0.7102 3.3248 4.9023
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.1817 4.6761 5.171 0.00129 **
## trend -2.1060 1.9529 -1.078 0.31662
## I(trend^2) 0.9216 0.1730 5.326 0.00109 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.976 on 7 degrees of freedom
## Multiple R-squared: 0.9812, Adjusted R-squared: 0.9758
## F-statistic: 182.5 on 2 and 7 DF, p-value: 9.137e-07cholesterol_data_forecast_table <- forecast(cholesterol_lm)
cholesterol_data_forecast <- cholesterol_data_forecast_table$fitted
accuracy(cholesterol_data_forecast, cholesterol_data)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 0 3.326378 2.872576 -0.4662172 6.287116 -0.3101741 0.3657351Finally we will plot the actual data and the forecasted data together.
cholesterol_data_forecast_m <- cbind(cholesterol_data, cholesterol_data_forecast)
cholesterol_data_forecast_m <- na.omit(cholesterol_data_forecast_m)
cholesterol_data_forecast_m## Time Series:
## Start = 1
## End = 10
## Frequency = 1
## cholesterol_data cholesterol_data_forecast
## 1 23.1 22.99727
## 2 21.3 23.65606
## 3 27.4 26.15803
## 4 34.6 30.50318
## 5 33.8 36.69152
## 6 43.2 44.72303
## 7 59.5 54.59773
## 8 64.4 66.31561
## 9 74.2 79.87667
## 10 99.3 95.28091plot.ts(cholesterol_data_forecast_m , type = "o", plot.type = "single", xlab = "Year", ylab = "Cholesterol", col = c("red", "blue"), main = "Quadratic Trend Forecast", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))
3. Seasonal Pattern
1. Seasonality Without Trend
Now we will start looking into Seasonal Pattern’s data. Let’s load a seasonal data from the Anderson book. As usual, first we will load the data, create a time series object and plot it
raw_data <- read_excel("/home/tanvir/Documents/ownCloud/Git_Repos/EWU_repos/3_Fall_2023/eco_204/ewu-eco204.github.io/pdf_files/slides/06_time_series/Codes/Umbrella.xlsx")
# Note that we used frequency = 4 because the data is quarterly data.
umbrella_data <- ts(raw_data$Sales, frequency = 4, start = c(2018,1))
plot.ts(umbrella_data, type = "o", xlab = "Year", ylab = "Sales")
Now let’s fit a linear model with seasonal dummies. The function
seasonaldummy() from the package forecast can
be used to create seasonal dummies.
seasonaldummy(umbrella_data)## Q1 Q2 Q3
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1
## [4,] 0 0 0
## [5,] 1 0 0
## [6,] 0 1 0
## [7,] 0 0 1
## [8,] 0 0 0
## [9,] 1 0 0
## [10,] 0 1 0
## [11,] 0 0 1
## [12,] 0 0 0
## [13,] 1 0 0
## [14,] 0 1 0
## [15,] 0 0 1
## [16,] 0 0 0
## [17,] 1 0 0
## [18,] 0 1 0
## [19,] 0 0 1
## [20,] 0 0 0Finally we fit a linear model with seasonal dummies.
umbrella_lm <- tslm(umbrella_data ~ season)
summary(umbrella_lm)##
## Call:
## tslm(formula = umbrella_data ~ season)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.0 -8.0 1.0 7.5 14.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 124.000 5.065 24.484 4.15e-14 ***
## season2 28.000 7.162 3.909 0.001249 **
## season3 -3.000 7.162 -0.419 0.680892
## season4 -29.000 7.162 -4.049 0.000931 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.32 on 16 degrees of freedom
## Multiple R-squared: 0.7989, Adjusted R-squared: 0.7611
## F-statistic: 21.18 on 3 and 16 DF, p-value: 8.104e-06Now we can calculate the forecast using the forecast()
function. Let’s see how we can do that.
Umrella_data_forecast_table <- forecast(umbrella_lm)
# we only need the point forecast now
Umrella_data_forecast <- Umrella_data_forecast_table$fitted
Umrella_data_forecast## Qtr1 Qtr2 Qtr3 Qtr4
## 2018 124 152 121 95
## 2019 124 152 121 95
## 2020 124 152 121 95
## 2021 124 152 121 95
## 2022 124 152 121 95accuracy(Umrella_data_forecast, umbrella_data)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -2.842225e-15 10.12917 8.9 -0.7581649 7.570416 0.5589669 0.3616559Finally we will plot the actual data and the forecasted data together.
# Nothing to omit here
Umrella_data_forecast_m <- cbind(umbrella_data, Umrella_data_forecast)
plot.ts(Umrella_data_forecast_m , type = "o", plot.type = "single", xlab = "Year", ylab = "Sales", col = c("red", "blue"), main = "Seasonal Trend Forecast", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))
1. Seasonality With Trend
We can also do seasonality and trend together. Here we will use the Smartphones data from the Anderson book.
raw_data <- read_excel("/home/tanvir/Documents/ownCloud/Git_Repos/EWU_repos/3_Fall_2023/eco_204/ewu-eco204.github.io/pdf_files/slides/06_time_series/Codes/SmartphoneSales.xlsx")
smartphone_data <- ts(raw_data$`Sales (1000s)`, frequency = 4, start = c(1,1))
smartphone_data## Qtr1 Qtr2 Qtr3 Qtr4
## 1 4.8 4.1 6.0 6.5
## 2 5.8 5.2 6.8 7.4
## 3 6.0 5.6 7.5 7.8
## 4 6.3 5.9 8.0 8.4plot.ts(smartphone_data, type = "o", xlab = "Year", ylab = "Sales")
Now we will fit a linear model with seasonal dummies and trend. Let’s see how we can do that.
smartphone_lm <- tslm(smartphone_data ~ trend + season)
summary(smartphone_lm)##
## Call:
## tslm(formula = smartphone_data ~ trend + season)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2988 -0.1569 -0.0075 0.1150 0.3663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.70563 0.13756 34.207 1.60e-12 ***
## trend 0.14562 0.01211 12.023 1.14e-07 ***
## season2 -0.67063 0.15368 -4.364 0.00113 **
## season3 1.05875 0.15511 6.826 2.85e-05 ***
## season4 1.36312 0.15745 8.657 3.06e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2167 on 11 degrees of freedom
## Multiple R-squared: 0.9763, Adjusted R-squared: 0.9676
## F-statistic: 113.2 on 4 and 11 DF, p-value: 7.376e-09Now we will calculate the forecast using the forecast()
function. Let’s see how we can do that.
smartphone_data_forecast_table <- forecast(smartphone_lm)
# we only need the point forecast now
smartphone_data_forecast <- smartphone_data_forecast_table$fitted
smartphone_data_forecast## Qtr1 Qtr2 Qtr3 Qtr4
## 1 4.85125 4.32625 6.20125 6.65125
## 2 5.43375 4.90875 6.78375 7.23375
## 3 6.01625 5.49125 7.36625 7.81625
## 4 6.59875 6.07375 7.94875 8.39875accuracy(smartphone_data_forecast, smartphone_data)## ME RMSE MAE MPE MAPE ACF1
## Test set 5.551115e-17 0.1796481 0.141875 -0.1032015 2.450953 0.3913853
## Theil's U
## Test set 0.1458527We can also plot the actual data and the forecasted data together.
# Nothing to omit here
smartphone_data_forecast_m <- cbind(smartphone_data, smartphone_data_forecast)
plot.ts(smartphone_data_forecast_m , type = "o", plot.type = "single", xlab = "Year", ylab = "Sales", col = c("red", "blue"), main = "Seasonal Trend Forecast", lwd = c(2, 1))
legend("topleft", legend=c("Actual", "Forecast"),col=c("red", "blue"), cex=0.8, lwd = c(2, 1))