#####################
#Time Series Analysis

######################################################################
#Time Domain

#Forecast library
library(forecast)

#White Noise
N=1000
mu=0
sigma=1

set.seed(5301)
z=rnorm(N,mu,sigma)

plot(1:N,z,type="l",xlab="t")


#Moving Average Process MA(q)

z.ma=rnorm(N,mu,sigma)

x=rep(0,N)

for(t in 3:N){
  x[t]=1*z.ma[t]+0.8*z.ma[t-1]+0.6*z.ma[t-2]
}

x.ma=x
plot(1:N,
     x.ma,
     type="l",
     xlab="t",
     ylab="x",
     main="MA(2) Process")





#Autocorrelations for z and x.ma

acf(z)

acf(x.ma,main="MA(2) Process")



#Random Walk

z.rw=2*rbinom(N,size=1,prob=.5)-1
z.rw

x=rep(0,N)

for(t in 2:N){
  x[t]=x[t-1]+z.rw[t]
}

x.rw=x
plot(1:N,
     x.rw,
     type="l",
     xlab="t",
     ylab="x",
     main="Random Walk AR(1)")


#An AR(2) Model

z.ar=rnorm(N,mu,sigma)
x=rep(0,N)

for(t in 3:N){
  x[t]=0.6*x[t-1]+0.3*x[t-2]+1*z.ar[t]
}

x.ar=x

plot(1:N,
     x.ar,
     type="l",
     xlab="t",
     ylab="x",
     main="AR(2) Process")


#Autocorrelations for z, x.ma, and x.ar
acf(z)

acf(x.ma,main="MA(2) Process")

acf(x.ar,main="AR(2) Process")

acf(x.rw,main="Random Walk AR(1) Process")

#An ARMA(1,1) Process

z.arma=rnorm(N,mu,sigma)
x=rep(0,N)

for(t in 2:N){
  x[t]=0.6*x[t-1]+1*z.arma[t]+0.4*z.arma[t-1]
}

x.arma=x

plot(1:N,
     x.arma,
     type="l",
     xlab="t",
     ylab="x",
     main="ARMA(1,1) Process")



######################################################
#Fitting ARIMA(p,d,q) Models

#MA(2)
#x[t]=1*z[t]+0.8*z[t-1]+0.6*z[t-2]
model.ma=auto.arima(x.ma)
summary(model.ma)

#AR(2)
#x[t]=0.6*x[t-1]+0.3*x[t-2]+1*z[t]
model.ar=auto.arima(x.ar)
summary(model.ar)


#ARMA(1,1)
#x[t]=0.6*x[t-1]+1*z[t] + 0.4*z[t-1
model.arma=auto.arima(x.arma)
summary(model.arma)


#Random Walk AR(1)
#x[t]=x[t-1]+z[t]
model.rw=auto.arima(x.rw)
summary(model.rw)




###################################################################
#Forecasting for ARIMA Models
###################################################################



###################################################################
#MA(2) Forecast

#Lead time
h=500

#Forecast Function
ma.forecast=forecast(model.ma,h)
names(ma.forecast)
head(ma.forecast$mean)
head(ma.forecast$upper)
head(ma.forecast$lower)


#Plotting the Forecast
plot(1:N,
     x.ma,
     xlim=c(1,N+h),
     type="l",
     xlab="t",
     ylab="x",
     main="MA(2) Process")

lines((N+1):(N+h),
      ma.forecast$mean,
      col="blue")

lines((N+1):(N+h),
      ma.forecast$upper[,2],
      col="blue")

lines((N+1):(N+h),
      ma.forecast$lower[,2],
      col="blue")

#Monte Carlo Simulations

ma.plot.wrapper=function(){
  plot(1:N,
       x.ma,
       xlim=c(1,N+h),
       type="l",
       xlab="t",
       ylab="x",
       main="MA(2) Process")
  
  lines((N+1):(N+h),
        ma.forecast$mean,
        col="blue")
  
  lines((N+1):(N+h),
        ma.forecast$upper[,2],
        col="blue")
  
  lines((N+1):(N+h),
        ma.forecast$lower[,2],
        col="blue")
}
ma.plot.wrapper()

ma.MC.wrapper=function(){
  x.fc=c(x.ma,rep(0,h))
  z.fc=c(z.ma,rnorm(h,mu,sigma))
  
  for(t in (N+1):(N+h)){
    x.fc[t]=1*z.fc[t]+0.8*z.fc[t-1]+0.6*z.fc[t-2]
  }
  
  lines((N+1):(N+h),
        x.fc[(N+1):(N+h)],
        col="red")
  
}
ma.MC.wrapper()


ma.plot.wrapper()
ma.MC.wrapper()



###################################################################
#AR(2) Forecast

#Lead time
h=500

#Forecast Function
ar.forecast=forecast(model.ar,h)

ar.plot.wrapper=function(){
  plot(1:N,
       x.ar,
       xlim=c(1,N+h),
       type="l",
       xlab="t",
       ylab="x",
       main="AR(2) Process")
  
  lines((N+1):(N+h),
        ar.forecast$mean,
        col="blue")
  
  lines((N+1):(N+h),
        ar.forecast$upper[,2],
        col="blue")
  
  lines((N+1):(N+h),
        ar.forecast$lower[,2],
        col="blue")
}
ar.plot.wrapper()


#Monte Carlo Simulations
ar.MC.wrapper=function(){
  x.fc=c(x.ar,rep(0,h))
  z.fc=c(z.ar,rnorm(h,mu,sigma))
  
  for(t in (N+1):(N+h)){
    x.fc[t]=0.6*x.fc[t-1]+0.3*x.fc[t-2]+1*z.fc[t]
  }
  
  lines((N+1):(N+h),
        x.fc[(N+1):(N+h)],
        col="red")
  
}
ar.MC.wrapper()


ar.plot.wrapper()
ar.MC.wrapper()



###################################################################
#Random Walk AR(1) Forecast

#Lead time
h=500

#Forecast Function
rw.forecast=forecast(model.rw,h)

rw.plot.wrapper=function(){
  upper.bound=max(c(x.rw,rw.forecast$upper[,2]))
  lower.bound=min(c(x.rw,rw.forecast$lower[,2]))
  
  plot(1:N,
       x.rw,
       xlim=c(1,N+h),
       ylim=c(lower.bound,upper.bound),
       type="l",
       xlab="t",
       ylab="x",
       main="Random Walk AR(1) Process")
  
  lines((N+1):(N+h),
        rw.forecast$mean,
        col="blue")
  
  lines((N+1):(N+h),
        rw.forecast$upper[,2],
        col="blue")
  
  lines((N+1):(N+h),
        rw.forecast$lower[,2],
        col="blue")
}
rw.plot.wrapper()


#Monte Carlo Simulations
rw.MC.wrapper=function(){
  x.fc=c(x.rw,rep(0,h))
  z.fc=c(z.rw,2*rbinom(h,size=1,prob=.5)-1)
  
  for(t in (N+1):(N+h)){
    x.fc[t]=x.fc[t-1]+z.fc[t]
  }
  
  lines((N+1):(N+h),
        x.fc[(N+1):(N+h)],
        col="red")
  
}
rw.MC.wrapper()

rw.plot.wrapper()
rw.MC.wrapper()


#######################################################################
#Dow Jones Data
#
#1/29/1985 to 11/6/2020

head(dowjones)

x=dowjones$Close
N=length(x)

plot(1:N,
     x,
     type="l",
     xlab="t",
     main="Dow Jones Industrial Average")


#Logarithmic Transformation
x=log(dowjones$Close)

plot(1:N,
     x,
     type="l",
     xlab="t",
     ylab="log(x)",
     main="Dow Jones Industrial Average")


N=ceiling(2/3*N)
h=length(x)-N
N
h


#Fitting an ARIMA Model

dow.model=auto.arima(x[1:N])
summary(dow.model)


#Forecast
dow.forecast=forecast(dow.model,h)

dow.plot.wrapper=function(){
  x.dow=x[1:N]
  upper.bound=max(c(x.dow,dow.forecast$upper[,2]))
  lower.bound=min(c(x.dow,dow.forecast$lower[,2]))
  
  plot(1:N,
       x.dow,
       xlim=c(1,N+h),
       ylim=c(lower.bound,upper.bound),
       type="l",
       xlab="t",
       ylab="log(x)",
       main="Dow Jones Industrial Average")
  
  lines((N+1):(N+h),
        dow.forecast$mean,
        col="blue")
  
  lines((N+1):(N+h),
        dow.forecast$upper[,2],
        col="blue")
  
  lines((N+1):(N+h),
        dow.forecast$lower[,2],
        col="blue")
}

dow.plot.wrapper()
lines((N+1):(N+h),
      x[(N+1):(N+h)],
      col="red")


#Removing Logarithmic Transformation

exp.dow.plot.wrapper=function(){
  x.dow=x[1:N]
  upper.bound=exp(max(c(x.dow,dow.forecast$upper[,2])))
  lower.bound=exp(min(c(x.dow,dow.forecast$lower[,2])))
  
  plot(1:N,
       exp(x.dow),
       xlim=c(1,N+h),
       ylim=c(lower.bound,upper.bound),
       type="l",
       xlab="t",
       ylab="x",
       main="Dow Jones Industrial Average")
  
  lines((N+1):(N+h),
        exp(dow.forecast$mean),
        col="blue")
  
  lines((N+1):(N+h),
        exp(dow.forecast$upper[,2]),
        col="blue")
  
  lines((N+1):(N+h),
        exp(dow.forecast$lower[,2]),
        col="blue")
  
  lines((N+1):(N+h),
        exp(x[(N+1):(N+h)]),
        col="red")
}
exp.dow.plot.wrapper()

#Quick Check
lines(1:nrow(dowjones),dowjones$Close,col="green")



#########################################################################
#Frequency Domain


library(freqdom)


#Deterministic Periodic Time Series
N=30
t=1:N

#Time Series Parameters
w=pi/3
theta=pi/6
R=2

x=R*cos(w*t+theta)
plot(t,x,type="l",ylim=c(-1.2*R,1.2*R))

t.range=seq(from=1,to=N,length=1000)
lines(t.range,R*cos(w*t.range+theta),col="red")


#Spectral Density
x.spec=spectral.density(x)

plot(x.spec$freq,x.spec$operators,type="l")
lines(c(w,w),c(0,3),col="blue")


#Restricting domain to [0,pi]
freq.range=seq(from=0,to=pi,length=1000)
x.spec=spectral.density(x,freq=freq.range)

plot(x.spec$freq,x.spec$operators,type="l")
lines(c(w,w),c(0,3),col="blue")



##################