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


#####  Reading Shiller & PPP data sets and defining variables

Sh_da <- read.csv("C://Financial Econometrics/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])
	

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


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


### Automatic Identification stock Returns: IBM, DIS, CAT
FX_da <- read.csv("http://www.bauer.uh.edu/rsusmel/4397/Stocks_FX_1973.csv",head=TRUE,sep=",")

x_date <- SFX_da$Date
x_cat <- SFX_da$CAT
x_dis <- SFX_da$DIS
x_ibm <- SFX_da$IBM
x_sp500 <- SFX_da$SP500
x_Mkt_RF <- SFX_da$Mkt_RF
x_SMB <- SFX_da$SMB
x_HML <- SFX_da$HML
x_RF <- SFX_da$RF

T <- length(x_cat)
T
lr_cat <- log(x_cat[-1]/x_cat[-T])
lr_dis <- log(x_dis[-1]/x_dis[-T])
lr_ibm <- log(x_ibm[-1]/x_ibm[-T])
lr_sp <- log(x_sp500[-1]/x_sp500[-T])
date_1 <- x_date[-1]
Mkt_RF <- x_Mkt_RF[-1]/100
SMB <- x_SMB[-1]/100
HML <- x_HML[-1]/100
RF <- x_RF[-1]/100

ibm_x <- lr_ibm - RF
cat_x <- lr_cat - RF
dis_x <- lr_dis - RF

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

acf(ibm_x)
pacf(ibm_x)

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

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


### Automatic Identification: Gold Price Changes, Oil Price Changes, US Inflation and Changes in FX rates

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

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

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


### Estimation & Diagnostic Testing: 
fit_arima_e_p10 <- arima(e_usd, order=c(0,0,10))
fit_arima_e_p10
autoplot(fit_arima_e_p10)			# Check roots (actually, inverse roots)
checkresiduals(fit_arima_e_p10)			# Check residuals (should look like White Noise, uncorrelated)	


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

fit_arima_I_p5 <- arima(us_I, order=c(1,1,5))
fit_arima_I_p5
autoplot(fit_arima_I_p5)
checkresiduals(fit_arima_I_p5)

fit_arima_I_p12 <- arima(us_I, order=c(1,1,12))
fit_arima_I_p12
autoplot(fit_arima_I_p12)
checkresiduals(fit_arima_I_p12)

fit_arima_I_p512 <- arima(us_I, order=c(5,1,12))
fit_arima_I_p512
autoplot(fit_arima_I_p512)
checkresiduals(fit_arima_I_p512)