############# Time Series: SES, HW ES & Forecasting (Code #18) #############

Car_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/TOTALNSA.csv", head=TRUE, sep=",")
names(Car_da)

x_car <- Car_da$TOTALNSA

### Stabilizing Transformations
library(tseries)
ts_car <- ts(x_car,start=c(1976,1),frequency=12) 
plot.ts(ts_car, xlab="Time", ylab="div", main="Total U.S. Vehicle Sales")
mean(x_car)
sd(x_car)

## Log Transformation
l_car <- log(ts_car)
plot.ts(l_car, xlab="Time", ylab="div", main="Total U.S. Vehicle Sales")
mean(l_car)
sd(l_car)

### SES Log Car Sales
mod_car <- HoltWinters(l_car, beta=FALSE, gamma=FALSE)
mod_car


### SES Log Dividend Changes
Sh_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/Shiller_data.csv", head=TRUE, sep=",")
names(Sh_da)
x_P <- Sh_da$P
x_D <- Sh_da$D
x_i <- Sh_da$Long_i
T <- length(x_P)
lr_p <- log(x_P[-1]/x_P[-T])
lr_d <- log(x_D[-1]/x_D[-T])
mod1 <- HoltWinters(lr_d[1:1750], beta=FALSE, gamma=FALSE)
mod1

plot(mod1$fitted, main= "Log Dividend Changes: SES" )

## One step forecast errors
T_last <- nrow(mod1$fitted)		# number of in-sample forecasts
h <- 25					# forecast horizon
ses_f <- matrix(0,h,1)			# Vector to collect forecasts
alpha <- 0.29
y <- lr_d			
T <- length(lr_d)
sm <- matrix(0,T,1)
T1 <- T - h + 1				# Start of forecasts
a <- T1					# index for while loop
sm[a-1] <- mod1$fitted[T_last]		# last in-sample forecast	
while (a <= T) {
  sm[a] = alpha * y[a-1] +  (1-alpha) * sm[a-1]
  a <- a + 1
} 

ses_f <- sm[T1:T]
ses_f
f_error_ses <- sm[T1:T] - y[T1:T]	# forecast errors
MSE_ses <- sum(f_error_ses^2)/h		# MSE
plot(ses_f, type="l", main ="SES Forecasts: Changes in Dividends")


forecast(mod1, h=25, level=.95)



### Holt-Winters
components_car <- decompose(l_car)
plot(components_car)

HW_car <- HoltWinters(l_car, seasonal = "additive")
HW_car

# Custom Holt-Winters fitting (Specifying alpha, beta & gamma)
HW_car_2 <- HoltWinters(l_car, alpha=0.2, beta=0.1, gamma=0.1)
HW_car_2

#Visual evaluation of fit
plot(l_car, main = "Log Car Sales", ylab="Log Car Sales")
lines(HW_car$fitted[,1], lty=2, col="blue")
lines(HW_car_2$fitted[,1], lty=2, col="red")

# Forecasting
HW_car.pred <- predict(HW_car, 24, prediction.interval = TRUE, level=0.95)

forecast(HW_car, h=24, level=.95)

# Visual evaluation of Fit and Forecasts
plot(l_car, main = "Forecasting Log Car Sales", ylab="Log CAr Sales")
lines(HW_car$fitted[,1], lty=2, col="blue")		# History of Series
lines(HW_car.pred[,1], col="red")			# Predictions
lines(HW_car.pred[,2], lty=2, col="orange")		# Upper 95% C.I. band
lines(HW_car.pred[,3], lty=2, col="orange") 		# Lower 95% C.I. band

# In-sample MSE
HW_car_mse <- sum(HW_car$fitted[,1]^2)/length(HW_car$fitted[,1])
HW_car_mse

# Multiplicative Model and Forecast
HW_car_m <- HoltWinters(l_car, seasonal ="multiplicative")
HW_car_m

HW_car_m.pred <- predict(HW_car_m, 24, prediction.interval = TRUE, level=0.95)

forecast(HW_car_m, h=24, level=.95)

#Visually evaluate the prediction
plot(l_car, main = "Forecasting Log Car Sales", ylab="Log Car Sales")
lines(HW_car_m$fitted[,1], lty=2, col="blue")		# History of Series
lines(HW_car_m.pred[,1], col="red")			# Predictions
lines(HW_car_m.pred[,2], lty=2, col="orange")		# Upper 95% C.I. band
lines(HW_car_m.pred[,3], lty=2, col="orange") 		# Lower 95% C.I. band

# In-sample MSE
HW_car_m_mse <- sum(HW_car_m$fitted[,1]^2)/length(HW_car_m$fitted[,1])
HW_car_m_mse


# Using R package forecast
library(forecast)
HW_car_f <- forecast(HW_car, h=24, level=c(80,95))

# Visualize forcast:
plot(HW_car_f, main = "Forecasting Log Car Sales", ylab="Log Car Sales") 
lines(HW_car_f$fitted, lty=2, col="purple")

# Evaluation of HW Model
HW_car_fe <- HW_car_f$residuals			#In-sample prediction errors

acf(HW_car_fe, lag.max=20, na.action=na.pass)

Box.test(HW_car_fe, lag=20, type="Ljung-Box")

hist(HW_car_fe)

# Testing accuracy of forecasts (Very informal)
dm.test(HW_car.pred[,1], HW_car_m.pred[,1], power=2)

dm.test(HW_car.pred[,1], HW_car_m.pred[,1], power=1)