############# Time Series - ARMA Processes: Identification & Estimation (Code #15 - Lecture 21) #############

##### Time Series - ARMA Processes: Identification of Simulated Data

### Simulation of MA process. Plots of ACF

## MA with Theta = 0.5
sim_ma1_5 <- arima.sim(list(order=c(0,0,1), ma=0.5), n=200) 
plot(arima.sim(list(order=c(0,0,1), ma=0.5), n=200), ylab="ACF",
main=(expression(MA(1)~~~theta==+.5)))

acf_ma1_5 <- acf(sim_ma1_5, main=(expression(MA(1)~~~theta==+.5)))
acf_ma1_5

T <- length(sim_ma1_5)
SE_ma <- sqrt(1/T)
SE_ma

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

## MA with Theta = -0.9
sim_ma1_9 <- arima.sim(list(order=c(0,0,1), ma=-0.9), n=200) 
plot(sim_ma1_9, ylab="ACF", main=(expression(MA(1)~~~theta==-.9)))

acf_ma1_9 <- acf(sim_ma1_9, main=(expression(MA(1)~~~theta==-.9)))
acf_ma1_9

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

## MA with Theta = 2
sim_ma1_2 <- arima.sim(list(order=c(0,0,1), ma=-2), n=200) 
acf_ma1_2 <- acf(sim_ma1_2, main=(expression(MA(1)~~~theta==-2)))
acf_ma1_2

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


### Simulation of AR process. Plots of ACF
## AR with Phi = 0.5
sim_ar1_5 <- arima.sim(list(order=c(1,0,0), ar=0.5), n=200) 
plot(sim_ar1_5, ylab="ACF", main=(expression(AR(1)~~~phi==+.5)))

acf_ar1_5 <- acf(sim_ar1_5, main=(expression(AR(1)~~~phi==+.5)))
acf_ar1_5

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

## AR with Phi = -0.9
sim_ar1_9 <- arima.sim(list(order=c(1,0,0), ar=-0.9), n=200) 
plot(sim_ma1_9, ylab="ACF", main=(expression(AR(1)~~~phi==-.9)))

acf_ar1_9 <- acf(sim_ar1_9, main=(expression(AR(1)~~~phi==-.9)))
acf_ar1_9

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

## AR with Phi = 2 (Can't be simulated by R, it's not stationary)
sim_ar1_2 <- arima.sim(list(order=c(1,0,0), ar=2), n=200) 

T_sim <- 200
u <- rnorm(T_sim, sd=1)				# Draw T_sim normally distributed errors
y_sim <- matrix(0,T_sim,1)			# vector to accumulate simulated data		
phi1 <- 2					# Change to create different correlation patterns
a <- k+1					# Time index for observations
mu <- 0.1
while (a <= T_sim) {
  y_sim[a] = mu + phi1 * y_sim[a-1] + u[a]	# y_sim simulated autocorrelated values	 
  a <- a + 1
} 
  y_sim_ar1 <- y_sim[(k+1):T_sim]

plot(y_sim_ar1, type="l", main =(expression(AR(1)~~~phi==+2)), ylab = "Simulated Series")
acf(y_sim_ar1)

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


## PACF - Simulated AR(2) 
sim_ar22 <- arima.sim(list(order=c(2,0,0), ar=c(0.5, 0.3)), n=200) 	#simulate AR(2) series
plot(sim_ar22, ylab="Simulated Series", main=(expression(AR(2))))
pacf_ar22 <- pacf(sim_ar22) 	

## PACF - Simulated MA(1) 
sim_ma1 <- arima.sim(list(order=c(0,0,1), ma = 0.5), n=200)  	
pacf(sim_ma1, main=(expression(MA(1)~~~theta==0.5)))

## PACF - Simulated ARMA(1,1) 
sim_arma11 <- arima.sim(list(order=c(1,0,1), ar=0.4, ma=0.5), n=200)  	
pacf(sim_arma11, main=(expression(ARMA(1,1)~~~phi==0.4~~~theta==0.5)))



##### Time Series - ARMA Processes: Identification of Real 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])

acf_p <- acf(lr_p)
pacf_p <- pacf(lr_p)

acf_p <- acf(lr_d)
pacf_p <- pacf(lr_d)

## Automatic Identification using Minic	(Two packages: forecast & caschrono)
library(caschrono)			# Install package caschrono (May not work for current version of R
armaselect(lr_p)

library(forecast)			# Install package forecast
auto.arima(lr_p)
auto.arima(lr_p, ic="bic", trace=TRUE)

auto.arima(lr_d)
auto.arima(lr_d, ic="bic", trace=TRUE)

library(weakARMA)			# Install package weakARMA (it gives you a lot of IC to choose from)
ar_p <- ARMA.selec(lr_p,P=4,Q=4)	# It may take a while with P & Q > 3
ar_p

ar_d <- ARMA.selec(lr_d,P=4,Q=4)
ar_d


##### Time Series - ARMA Processes: Identification of Real Data

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

x_date <- FMX_da$Date
us_CPI <- FMX_da$CPI_US
nor_CPI <- FMX_da$CPI_Norway
S_nok <- FMX_da$NOK_USD
S_gbp <- FMX_da$USD_GBP
x_Mkt <- FMX_da$Mkt_RF
x_SMB <- FMX_da$SMB
x_HML <- FMX_da$HML
x_RF <- FMX_da$RF
x_oil <- FMX_da$Crude_WTI_Oil
x_gold <- FMX_da$Gold
x_cons <- FMX_da$Cons_sent

T <- length(us_CPI)
us_I <- log(us_CPI[-1]/us_CPI[-T])
nor_I <- log(nor_CPI[-1]/nor_CPI[-T])
e_nok <- log(S_nok[-1]/S_nok[-T])
e_gbp <- log(S_gbp[-1]/S_gbp[-T])
oil <- log(x_oil[-1]/x_oil[-T])
gold <- log(x_gold[-1]/x_gold[-T])
cons <- log(x_cons[-1]/x_cons[-T])

acf_e <- acf(e_gbp, main="USD/GBP FX Rate: Monthly Log Changes (1973-205)")
acf_e

plot(e_gbp, main="USD/GBP FX Rate: Monthly Log Changes (1973-205)", type="l", col="blue")

pacf_e <- pacf(e_gbp)
pacf_e

acf_e_n <- acf(e_nok, main="NOK/USD FX Rate: Monthly Log Changes (1973-205)")
acf_e_n

pacf_e_n <- pacf(e_nok)
pacf_e_n

acf_inf <- acf(us_I)
acf_inf

pacf_inf <- pacf(us_I)
pacf_inf

acf_o <- acf(oil)
acf_o

pacf_o <- pacf(oil)
pacf_o

acf_g <- acf(gold)
acf_g

pacf_g <- pacf(gold)
pacf_g

acf_c <- acf(cons)
acf_c

pacf_c <- pacf(cons)
pacf_c

ar_e <- ARMA.selec(e_gbp, P=3,Q=3)
ar_e

ar_o <- ARMA.selec(oil,P=3,Q=3)
ar_o

ar_g <- ARMA.selec(gold, P=3,Q=3)
ar_g


### Non-Stationary Time Series: Detrend or Differentiation?

## Levels
plot(x_P, type="l", col="blue", ylab ="S&P 500 Prices", xlab ="Time")
title("U.S. Stock Prices")

plot(x_P, type="l", col="blue", ylab ="S&P 500 Dividends", xlab ="Time")
title("U.S. Stock Dividends")