##############
# ARMA Model #
##############

# white noise
set.seed(1234)
n = 100
e = rnorm(n)
plot(e, type="o", main="white noise")
e[1:5]

# MA(1)
theta = 0.6
y = rep(0,n)
y[1] = e[1]
for (t in 2:n) {y[t]=e[t]+theta*e[t-1]}

par(mfrow=c(2,2)) 
plot(y, type="o", main="MA1")
acf(y)
plot(e, type="o", main="white noise")
acf(e)

arima(y, order = c(0,0,1))

# exercise 
theta1 = 0.6; theta2 = 0.4
y = rep(0,n)
y[1] = e[1]; y[2] = e[2]
for (t in 3:n) {y[t]=e[t]+theta1*e[t-1]+theta2*e[t-2]}
acf(y)
arima(y, order = c(0,0,2))
arima(e, order = c(0,0,1))


# AR(1)
y = rep(0,n)
y[1] = e[1]
for (t in 2:n) {y[t]=0.6*y[t-1]+e[t]}
y1 = y
for (t in 2:n) {y[t]=1*y[t-1]+e[t]}
y2 = y

par(mfrow=c(2,3)) 
plot(y1, type="o", main="stationary AR(1)")
acf(y1)
pacf(y1)
plot(y2, type="o", main="unit root")
acf(y2)
pacf(y2)

# Estimate AR(1)
arima(y1, order = c(1,0,0))
arima(y2, order = c(1,0,0))

# Check adequacy
Box.test(arima(y1, order = c(1,0,0))$res, lag = 1, type = "Ljung")
Box.test(arima(y2, order = c(1,0,0))$res, lag = 1, type = "Ljung")

# plot y against yhat
y1hat = y1 - arima(y1, order = c(1,0,0))$res
y2hat = y2 - arima(y2, order = c(1,0,0))$res

par(mfrow=c(1,2))
plot(y1, type="o", main="stationary AR(1)")
lines(y1hat, lty=2)
plot(y2, type="o", main="unit root")
lines(y2hat, lty=2)

# Forecasting

# library(forecast)
# forecast(Arima(y1,order=c(1,0,0),include.drift=F),h=3)
# forecast(Arima(y2,order=c(1,1,0),include.drift=F),h=3)


# unit root test
library(tseries)
adf.test(y1, alternative = c("s"), k = 1)
adf.test(y2, alternative = c("s"), k = 1)

####################
# smoothing method #
####################

# user-defined function
fma = function(data, k) {
n = length(data)
oneside = (k-1)/2
fs = rep(NA, n)
for (t in (1+oneside):(n-oneside)) {
fs[t] = mean(data[(t-oneside):(t+oneside)])
}
return(fs)
}

# Simulate smoothing
set.seed(100)
t = 1:40
signal = 0.5*t+0.05*t^2
noise = 10*rnorm(40)
y = signal + noise
plot(t, y, type="l",col="black",lwd=3, main="Smoothing")
lines(signal,lty=c(1),col="blue",lwd=3)
lines(fma(y,7),lty=c(1),col="red",lwd=3)
legend(0, 90, c("y", "signal", "7-average"), col = c("black", "blue", "red"),lty = c(1, 1, 1), merge = TRUE)


# compare  k = 3 and k = 7
plot(t, y, type="l",col="black",lwd=3, main="Comparing k=7 (red) to k=3 (blue)")
lines(fma(y,7),lty=c(1),col="red",lwd=3)
lines(fma(y,3),lty=c(1),col="blue",lwd=3)

sd(fma(y,7),na.rm=T)
sd(fma(y,3),na.rm=T)
length(na.omit(fma(y,7)))
length(na.omit(fma(y,3)))


# user-defined ESF
fes = function(data, alpha) {
n = length(data)
yes = rep(NA, n)
yes[2] = y[1]
for (t in 2:n) {
yes[t+1]=yes[t] + alpha*(y[t]-yes[t])
}
return(yes)
}

# compare ESF with alpha = 0.9 and alpha = 0.3
plot(t, y, type="l",col="black",lwd=3, main="Comparing ESF with alpha = 0.9 (red) to alpha = 0.3 (blue)")
lines(fes(y,0.9),lty=c(1),col="red",lwd=3)
lines(fes(y,0.3),lty=c(1),col="blue",lwd=3)


# covid-data
#setwd("C:/Users/lij14/Dropbox/400")
setwd("C:/Users/chnlj/Dropbox/400")
data = read.table("672_covidohio_data.txt", header=T)
case = data[,2]
par(mfrow=c(1,2)) 
plot(case, type="o", main="Whole sample")
plot(case[310:length(case)], type="o", main="Recent data")

# subsample 1/8/2021 to 3/7/2021
y = case[310:length(case)]
ts(y, frequency=365, start=c(2021, 1, 8))


# seven-day one-sided simple average (forecasting) and ESF
n = length(y)
yav = rep(NA, n)
for (t in 8:n) {yav[t]=sum(y[(t-7):(t-1)])/7}

plot(y, type="o", main="Smoothng forecasts for Recent Covid Cases in OH")
lines((1:n),yav,lty=2)
lines(fes(y,0.9), lty=3)
legend(40, 8000, c("y", "7-average", "ESF"), lty = c(1, 2, 3), merge = TRUE)

cat("7-average forecast for next period is ", yav[length(yav)], "\n")
cat("ESF for next period is ", fes(y,0.9)[length(fes(y,0.9))], "\n")