############# Time Series: Seasonalities & Forecasting (Code #17) #############

##### Seasonal Patterns: Monthly U.S. Retail Sales
AHB_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/Appl_Homes_Bread_Retail_2024.csv", head=TRUE, sep=",")
names(AHB_da)
summary(AHB_da)
x_A <- AHB_da$Appl_Man
x_H <- AHB_da$New_homes
x_B <- AHB_da$Bread_price
x_R <- AHB_da$Retail_Sales
T <- length(x_A)
lr_a <- log(x_A[-1]/x_A[-T])
lr_h <- log(x_H[-1]/x_H[-T])
lr_b <- log(x_B[-1]/x_B[-T])
lr_r <- log(x_R[-1]/x_R[-T])


z_lev <- x_R
lr_z <- lr_r


## Level
x.ts = ts(z_lev, frequency = 12, start=c(1992, 2))
plot.ts(x.ts, col="blue",ylab ="Retail Sales (USD)", main="Monthly Retail Sales (USD): 1992:Feb - 2024:Sep")

acf(z_lev)
pacf(lr_z)

T <- length(z_lev)
trend <- c(1:T)						# create trend
length(z_lev)
det_Z <- lm(z_lev ~ trend)				# regression to get detrended e
detrend_Z <- det_Z$residuals
plot(detrend_Z, type="l", col="blue", ylab ="Detrended Level", xlab ="Time")
title("Detrended Level")

# Add squared trend
trend2 <- trend^2
det_Z2 <- lm(z_lev ~ trend + trend2)		# regression to get detrended e
detrend_Z2 <- det_Z2$residuals
plot(detrend_Z2, type="l", col="blue", ylab ="Detrended Level with linear & quadratic trend", xlab ="Time")
title("Detrended Level")

# Work instead with log prices
l_Z <- log(z_lev)
det_lZ <- lm(l_Z ~ trend)			# regression to get detrended e
detrend_lZ <- det_lZ$residuals
plot(detrend_lA, type="l", col="blue", ylab ="Detrended Log Level", xlab ="Time")
title("Detrended Log Level")

det_lZ2 <- lm(l_Z ~ trend + trend2)		# regression to get detrended e
det_lZ2 <- det_lZ2$residuals
plot(det_lZ2, type="l", col="blue", ylab ="Det Log Manufactured Appliances", xlab ="Time")
title("Detrended Log Manufactured Appliances with linear and quadratic trends")

## Stochastic Trend: Differentiate
diff_Z <- diff(z_lev)
plot(diff_Z,type="l", col="blue", ylab ="Differenced Manufactured Appliances", xlab ="Time")
title("Differenced Manufactured Appliances")

## Try % change
plot(lr_z,type="l", col="blue", ylab ="Percentage Chage Manufactured Appliances", xlab ="Time")
title("% Chage Manufactured Appliances")

acf(lr_z)
pacf(lr_z)


###### Removing Seasonal Pattern
## Create seasonal dummies (Monthly effects)
zz <- lr_z
T <- length(lr_z)
M1 <- rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create January dummy
M2 <- rep(c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create March dummy
M3 <- rep(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create April dummy
M4 <- rep(c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create May dummy
M5 <- rep(c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create June dummy
M6 <- rep(c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create Jul dummy
M7 <- rep(c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create Aug dummy
M8 <- rep(c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0), (length(zz)/12+1))	# Create Sep dummy
M9 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), (length(zz)/12+1))	# Create Oct dummy
M10 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0), (length(zz)/12+1))	# Create Oct dummy
M11 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0), (length(zz)/12+1))	# Create Oct dummy
M12 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), (length(zz)/12+1))	# Create Oct dummy

Mar <- M1[1:T]
Apr <- M2[1:T]
May <- M3[1:T]
Jun <- M4[1:T]
Jul <- M5[1:T]
Aug <- M6[1:T]
Sep <- M7[1:T]
Oct <- M8[1:T]
Nov <- M9[1:T]
Dec <- M10[1:T]
Jan <- M11[1:T]

Seas <- cbind(Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec, Jan)
z_seas <- lm(lr_z ~ Seas)
summary(z_seas)
filt_z <- z_seas$residuals

plot(filt_z, type = "l")
acf(filt_z)
pacf(filt_z)

lr_z <- filt_z

## ARIMA Model: Auto.arima
auto.arima(lr_z, ic="bic", trace=TRUE)


auto_z <- auto.arima(lr_z, ic="bic", trace=TRUE)
autoplot(auto_z)
checkresiduals(auto_z)


#### Variance Stabilizing Transformation  

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)



## Real Estate Example: LA Home Prices with Seasonal Pattern
RE_da <- read.csv("https://www.bauer.uh.edu/rsusmel/4397/Real_Estate_2022.csv", head=TRUE, sep=",")
x_la <- RE_da$LA_c
zz <- x_la
T <- length(zz)
plot(x_la, type="l", main="Changes in Log Real Estate Prices in LA")

acf(x_la)
pacf(x_la)

Feb1 <- rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create January dummy
Mar1 <- rep(c(0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create March dummy
Apr1 <- rep(c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create April dummy
May1 <- rep(c(0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create May dummy
Jun1 <- rep(c(0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create June dummy
Jul1 <- rep(c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create Jul dummy
Aug1 <- rep(c(0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), (length(zz)/12+1))	# Create Aug dummy
Sep1 <- rep(c(0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0), (length(zz)/12+1))	# Create Sep dummy
Oct1 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0), (length(zz)/12+1))	# Create Oct dummy
Nov1 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0), (length(zz)/12+1))	# Create Oct dummy
Dec1 <- rep(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0), (length(zz)/12+1))	# Create Oct dummy
seas1 <- cbind(Feb1, Mar1, Apr1, May1, Jun1, Jul1, Aug1, Sep1, Oct1, Nov1, Dec1)
seas <- seas1[1:T,]

x_la_fit_sea <- lm(x_la ~ seas)				# Regress x_la against constant + seasonal dummies
summary(x_la_fit_sea)

x_la_filt <- x_la_fit_sea$residuals			# residuals, et = filtered x_la series

library(forecast)
fit_ar_la_filt <- auto.arima(x_la_filt)			# use auto.arima to look for a good model
fit_ar_la_filt

checkresiduals(fit_ar_la_filt)
autoplot(fit_ar_la_filt)



### Smoothing & Forecasting
library(tseries)
library(forecast)

## CAR SALES: 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
##  Reading Shiller Data

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)