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

##### 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
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)
	


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_OIL
x_gold <- FMX_da$GOLD

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])

acf_e <- acf(e_gbp)
acf_e

pacf_e <- pacf(e_gbp)
pacf_e

acf_inf <- acf(us_I)
acf_inf

pacf_inf <- pacf(us_I)
pacf_inf

acf_g <- acf(gold)
acf_g

pacf_g <- pacf(gold)
pacf_g


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

## Trend stationarity: Detrend
T <- length(x_P)					# length of series
trend <- c(1:T)						# create trend
det_P <- lm(x_P ~ trend)				# regression to get detrended e
detrend_P <- det_P$residuals
plot(detrend_P, type="l", col="blue", ylab ="Detrended U.S. Prices", xlab ="Time")
title("Detrended U.S. Stock Prices")

# Add squared trend
trend2 <- trend^2
det_P <- lm(x_P ~ trend + trend2)		# regression to get detrended e
detrend_P <- det_P$residuals
plot(detrend_P, type="l", col="blue", ylab ="Detrended U.S. Prices", xlab ="Time")
title("Detrended U.S. Stock Prices with linear and quadratic trends")

# Work instead with log prices
l_P <- log(x_P)
det_lP <- lm(l_P ~ trend)			# regression to get detrended e
detrend_lP <- det_lP$residuals
plot(detrend_lP, type="l", col="blue", ylab ="Detrended Log U.S. Prices", xlab ="Time")
title("Detrended Log U.S. Stock Prices")

det_lP2 <- lm(l_P ~ trend + trend2)		# regression to get detrended e
det_lP2 <- det_lP$residuals
plot(det_lP2, type="l", col="blue", ylab ="Det Log U.S. Prices", xlab ="Time")
title("Detrended Log U.S. Stock Prices with linear and quadratic trends")

## Stochastic Trend: Differentiate
diff_P <- diff(x_P)
plot(diff_P,type="l", col="blue", ylab ="Differenced U.S. Stock Prices", xlab ="Time")
title("Differenced U.S. Stock Prices")


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

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

## Trend stationarity: Detrend 
T <- length(x_D)					# length of series
trend <- c(1:T)						# create trend
det_D <- lm(x_D ~ trend)				# regression to get detrended e
detrend_D <- det_D$residuals
plot(detrend_D, type="l", col="blue", ylab ="Detrended U.S. Prices", xlab ="Time")
title("Detrended U.S. Stock Prices")

# Add squared trend
trend2 <- trend^2
det_D <- lm(x_D ~ trend + trend2)		# regression to get detrended e
detrend_D <- det_D$residuals
plot(detrend_D, type="l", col="blue", ylab ="Detrended U.S. Prices", xlab ="Time")
title("Detrended U.S. Stock Prices with linear and quadratic trends")

# Work instead with log dividends
l_D <- log(x_D)
det_lD <- lm(l_D ~ trend)			# regression to get detrended e
detrend_lD <- det_lD$residuals
plot(detrend_lD, type="l", col="blue", ylab ="Detrended Log U.S. Prices", xlab ="Time")
title("Detrended Log U.S. Stock Dividends")

det_lD2 <- lm(l_D ~ trend + trend2)		# regression to get detrended e
det_lD2 <- det_lD$residuals
plot(det_lD2, type="l", col="blue", ylab ="Det Log U.S. Dividends", xlab ="Time")
title("Detrended Log U.S. Stock Dividends with linear and quadratic trends")

## Stochastic Trend: Differentiate
diff_D <- diff(x_D)
plot(diff_D,type="l", col="blue", ylab ="Differenced U.S. Stock Dividends", xlab ="Time")
title("Differenced U.S. Stock Prices")

diff_det_D <- diff(detrend_D)
plot(diff_det_D,type="l", col="blue", ylab ="Differenced Detrended U.S. Stock Dividends", xlab ="Time")
title("Differenced U.S. Stock Dividends")

diff_det_lD <- diff(detrend_lD)
plot(diff_det_lD,type="l", col="blue", ylab ="Differenced Log Detrended U.S. Stock Dividends", xlab ="Time")
title("Differenced Log U.S. Stock Dividends")


## ARIMA? - Stock Prices, Stock Dividends & Interest Rates
acf_P <- acf(x_P)
acf_P

acf_D <- acf(x_D)
acf_D

acf_i <- acf(x_i)
acf_i


## Identification - Difference: Changes in Stock Prices & Changes in Stock Dividends
acf_p <- acf(lr_p)
pacf_p <- pacf(lr_p)

acf_d <- acf(lr_d)
pacf_d <- pacf(lr_d)


## Automatic Identification using Minic
library(caschrono)				# Need to install package caschrono first
armaselect(lr_p)

library(forecast)				# Need to install package forecast first
auto.arima(lr_p)
auto.arima(lr_p, ic="bic", trace=TRUE)