# load the data
rs <- read.csv("C:/Users/Sanjiv/Desktop/Math547/Week9/RadioSales.csv")
attach(rs)1
Sales <- ts(Sales, frequency = 4) # create a time series object
plot.ts(Sales, xlab="Year", main="Time Series Plot of C&D Marine Radios Sales")
We observe that sales increase from year to year. We also observe that 2nd quarter sales are usually the biggest, followed by 3rd quarter sales, while 4th quarter sales being the smallest, usually slightly smaller than 1st quarter sales.
2
# install.packages("forecast")
library(forecast)
ma <- ma(Sales, order=4, centre=F)
centered.ma <- ma(Sales, order=4, centre=T)
plot.ts(Sales, main="Moving Average", ylab="Sales", type="o")
lines(ma, col="blue", lwd=2, type="o")
lines(centered.ma, col="orange", lwd=2, type="o")
legend("topleft", lty=1, lwd=c(1,2,2), col=c(1,"blue","orange"),
legend=c("data","non-centered MA","centered MA"), pch=1)
# obtain seasonal indicies for each quarter
Sales/centered.ma ## Qtr1 Qtr2 Qtr3 Qtr4
## 1 NA NA 1.081081 0.395062
## 2 0.898876 1.484536 1.153846 0.482759
## 3 0.848485 1.434483 1.187097 0.592593
## 4 0.915663 1.287356 1.092896 0.750000
## 5 0.875622 1.314010 1.056604 0.777778
## 6 0.872727 1.303167 1.071429 0.689655
## 7 0.929461 1.264822 NA NAAs we can see, seasonal indices in order from high to low are: 2nd quarter, 3rd quarter, 1st quarter, 4th quarter. This makes sense since warmer seasons (2,3) are more appropriate for marine vacations than colder seasons (1,4).
3
# deseasonalize the time series
fit <- decompose(Sales, type="multiplicative") # decompose time series
deseason <- Sales/fit$seasonal # deseasonalized time series
deseason## Qtr1 Qtr2 Qtr3 Qtr4
## 1 6.67312 11.01581 8.94181 6.44279
## 2 11.12186 13.21898 13.41271 11.27488
## 3 15.57061 19.09408 20.56616 19.32836
## 4 21.13154 20.56285 22.35452 28.99254
## 5 24.46810 24.96918 25.03707 33.82463
## 6 26.69247 26.43795 26.82543 32.21394
## 7 31.14121 29.37550 31.29633 43.48881# computing deseasonalized trend line
trend.fit <- lm(deseason ~ Period)
trend.fit##
## Call:
## lm(formula = deseason ~ Period)
##
## Coefficients:
## (Intercept) Period
## 6.33 1.05fitted.deseason <- ts(fitted(trend.fit), frequency=4) # deseasonalized fitted values
fitted <- fitted.deseason*fit$seasonal
plot.ts(Sales, main="Time Series Decomposition", ylab="Sales", ylim=c(0,45))
lines(deseason, col="blue", lty=2, type="o")
lines(fitted.deseason, col="red")
lines(fitted, col="green", lty=2, type="o")
legend("topleft",lty=c(1,2,1,2), col=c(1, "blue","red","green"),
legend=c("original data", "deseasonalized data", "trend line", "fitted values"))
As we can see, there is increasing trend (represented by the red line with intercept 6.322 and slope 1.055).
4
# forecasts using decomposed time series
new.time <- data.frame(Period=c(29,30,31,32))
forecasts <- predict(trend.fit, new.time)*fit$figure
forecasts <- ts(forecasts, frequency=4, start=8)
forecasts## Qtr1 Qtr2 Qtr3 Qtr4
## 8 33.1928 51.7046 43.6443 24.8840plot.ts(Sales, main="Forecasts with Time Series Decomposition", ylab="Sales",
ylim=c(0,51), xlim=c(1,9))
lines(deseason, col="blue", lty=2, type="o")
lines(fitted.deseason, col="red")
lines(fitted, col="green", lty=2, type="o")
lines(forecasts, col="magenta", lty=2, type="o")
legend("topleft", lty=c(1,2,1,2,2), col=c(1, "blue","red","green", "magenta"),
legend=c("original data", "deseasonalized data", "trend line", "fitted values","forecasts"))
5
Here we will treat Quarter as factor, so essentially we fit four regression lines, one for each quarter.
lm.fit <- lm(Sales ~ Period + factor(Quarter), data = rs)
summary(lm.fit)##
## Call:
## lm(formula = Sales ~ Period + factor(Quarter), data = rs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.339 -1.016 0.085 0.730 3.170
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.9487 0.7879 6.28 2.1e-06 ***
## Period 0.9710 0.0382 25.41 < 2e-16 ***
## factor(Quarter)2 9.6004 0.8657 11.09 1.0e-10 ***
## factor(Quarter)3 4.2009 0.8682 4.84 7.0e-05 ***
## factor(Quarter)4 -4.9129 0.8724 -5.63 9.9e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.62 on 23 degrees of freedom
## Multiple R-squared: 0.975, Adjusted R-squared: 0.971
## F-statistic: 228 on 4 and 23 DF, p-value: <2e-16predict(lm.fit, newdata = data.frame(Period=c(29,30,31,32),Quarter=c(1,2,3,4)))## 1 2 3 4
## 33.1071 43.6786 39.2500 31.1071plot(Sales ~ Period, data=rs, type="l", main="Linear Regression Model")
abline(lm.fit$coefficients[1],lm.fit$coefficients[2], col=2, lty=2)
abline(lm.fit$coefficients[1]+lm.fit$coefficients[3],lm.fit$coefficients[2], col=3, lty=2)
abline(lm.fit$coefficients[1]+lm.fit$coefficients[4],lm.fit$coefficients[2], col=4, lty=2)
abline(lm.fit$coefficients[1]+lm.fit$coefficients[5],lm.fit$coefficients[2], col=5, lty=2)
legend("topleft",lty=2, col=c(2,3,4,5),
legend=c("Quarter 1","Quarter 2","Quarter 3","Quarter 4"))
6
# for time series model prediction is:
forecasts## Qtr1 Qtr2 Qtr3 Qtr4
## 8 33.1928 51.7046 43.6443 24.8840# for linear regression model prediction is:
predict(lm.fit, newdata = data.frame(Period=c(29,30,31,32),Quarter=c(1,2,3,4)))## 1 2 3 4
## 33.1071 43.6786 39.2500 31.1071We notice that predicted value for the 1st quarter (approx. 33.19281) is same for both models. For the 2nd and 3rd quarter time series model predicted higher sales than linear regression model, while for the 4th quarter regression predicted higher sales.
Calculating forecast accuracy on the historical data:
# for time series model
accuracy(fitted, Sales)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -0.231906 2.60001 1.98127 -2.63585 11.1275 0.159896 0.264677# for linear regression
accuracy(lm.fit$fitted.values, Sales)## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set -1.58566e-16 1.46635 1.12819 -0.746669 6.07603 0.425378 0.15561Mean absolute percentage error (MAPE) of 6.076028 is lower for the regression method.
7
# forecast error for time series model
35-33.19281## [1] 1.80719# forecast error for linear regression
35-33.10714## [1] 1.89286For the time series model we have smaller forecast error of 1.80719.