# Time Series Analysis Laboratory with R
# Prof. Lea Petrella
# Faculty of Economics
# Department of Methods and Models for Economics, Territory and Finance

# Sapienza University of Rome
# a.y. 2019-2020
#####################################

# Useful materials and links
# R: https://www.r-project.org
# RStudio: https://www.rstudio.com/home/
# Manual: The Art of R Programming
# Repositories & Code: GitHub (https://github.com)
#####################################


########################################
# Lesson 5: R Laboratory (Part 3)
########################################


##########################################################################
#
#   Index:
# 1)moving average processes
# 2)reversibility of a process
# 3)auto-regressive processes
# 4)concept of stationarity
# 5)the random walk process
#
################################################## 

**********************************
#####TIME SERIES SIMULATION WITH TREND AND ACF AND PACF VERIFICATION #####
# Choose parameters
alpha = 2.0
beta = 0.3

# Simulate a series of 200 observations
N = 200
Sim = c ()
for (i in 1: N) {
  Sim [i] = alpha + beta * i + rnorm (1, mean = 0, sd = 2)
}
par (mfrow = c (1,1))
plot.ts (Sim) # graphical analysis
abline (alpha, beta, col = "red") # the straight line "true"

# Study the autocorrelation function and the partial autocorrelation function
par (mfrow = c (2,1))
acf (Sim)
PACF (Sim)
# Does the ACF look like a White Noise process? How can we make it look like a White Noise?


HOMEWORK: estimate a regression model and analyse the residuals 
###############################

rm (list = ls ())
# The ARIMA process of class (p, i, q) can be simulated by using the 'arima.sim' command
? arima.sim

##### MOVING AVERAGE PROCESSES #####
# Let's start with the MA processes (1)
# simulate a process (invertible) with positive coefficient
n = 1000
Y1 = arima.sim (n = n, model = list (ma = c (0.8)), sd = sqrt (1), rand.gen = rnorm)
# and with negative coefficient
Y2 = arima.sim (n = n, model = list (ma = c (-0.8)), sd = sqrt (1), rand.gen = rnorm)

par (mfrow = c (2,1))
ts.plot (Y1, main = 'positive coefficient')
ts.plot (Y2, main = 'negative coefficient')

# compare the acf of the two processes
acf (Y1, main = 'positive coefficient')
acf (Y2, main = 'negative coefficient')

# compare the pacf of the two processes
pacf (Y1, main = 'positive coefficient')
pacf (Y2, main = 'negative coefficient')

# the same processes are obtained by simulating them without the command 'arima.sim'
# we initialize two vectors Ysim1 and Ysim2
Ysim1 = c ()
Ysim2 = c ()
# fix the first observation and generate the error component
Ysim1 [1] = 0
Ysim2 [1] = 0
shock = rnorm (n, 0, 1)

for (i in 2: n) {
Ysim1 [i] = 0.8 * shock [i-1] + shock [i]
Ysim2 [i] = -0.8 * shock [i-1] + shock [i]
}

par (mfrow = c (2,1))
ts.plot (Ysim1, main = 'positive coefficient')
ts.plot (Ysim2, main = 'negative coefficient')

# we compare the acf of the two processes
acf (Ysim1, main = 'positive coefficient')
acf (Ysim2, main = 'negative coefficient')

# we compare the pacf of the two processes
pacf (Ysim1, main = 'positive coefficient')
pacf (Ysim2, main = 'negative coefficient')

#Invertibility of a process:
# For every invertible MA (1) representation there is a non-invertible representation MA (1) with the same ACF.
theta = .8
theta.star = 1 / theta
Y3 = arima.sim (n = n, list (ma = theta.star), sd = sqrt (1))
Y4 = arima.sim (n = n, list (ma = -theta.star), sd = sqrt (1))


# we compare the acf of the two processes
acf (Y1, main = 'Reversible process') # reversible process
acf (Y3, main = 'Similar non-invertible process') # non-invertible process

# Compare the pacf of the two processes
pacf (Y1, main = 'Reversible process') # reversible process
pacf (Y3, main = 'Similar non-invertible process') # non-invertible process

# The two processes have identical properties and therefore it would be impossible to recognize them only via time series. This identification problem
# can be solved by constraining the parameter in the interval (-1, + 1).
# Such constraint also allows to rewrite the process as a stationary AR (1)

# Let's move on to the MA (2) processes ... there are four cases:
# 1. both positive coefficients
Y5 = arima.sim (n = n, list (ma = c (0.6, 0.3)), sd = sqrt (1))
# 2. both negative coefficients
Y6 = arima.sim (n = n, list (ma = c (-0.6, -0.3)), sd = sqrt (1))
# 3. coefficients with alternate marks
Y7 = arima.sim (n = n, list (ma = c (0.6, -0.3)), sd = sqrt (1))
# 4. coefficients with alternate marks
Y8 = arima.sim (n = n, list (ma = c (-0.6, 0.3)), sd = sqrt (1))

par (mfrow = c (2,2))
ts.plot (Y5, main = 'positive coefficients')
ts.plot (Y6, main = 'negative coefficients')
ts.plot (Y7, main = 'alternate coefficients (+, -)')
ts.plot (Y8, main = 'alternate coefficients (-, +)')

# Compare the acf of the four processes
acf (Y5, main = 'positive coefficients')
acf (Y6, main = 'negative coefficients')
acf (Y7, main = 'alternate coeff. (+, -)')
acf (Y8, main = 'alternate coeff. (-, +)')

# Compare the pacf of the four processes
pacf (Y5, main = 'positive coefficients')
pacf (Y6, main = 'negative coefficients')
pacf (Y7, main = 'alternate coeff. (+, -)')
pacf (Y8, main = 'alternate coeff. (-, +)')

# The same processes can be obtained by simulating them without the command 'arima.sim'
# we initialize the vectors that will contain the data
# NB! Now it is necessary to fix both the first and the second observation
Ysim5 = c (); Ysim5 [1] = 0; Ysim5 [2] = 0;
Ysim6 = c (); Ysim6 [1] = 0; Ysim6 [2] = 0;
Ysim7 = c (); Ysim7 [1] = 0; Ysim7 [2] = 0;
Ysim8 = c (); Ysim8 [1] = 0; Ysim8 [2] = 0;
for (i in 3:n){
  Ysim5[i] =   0.6 * shock[i-1] + 0.3 * shock[i-2] + shock[i]
  Ysim6[i] =  -0.6 * shock[i-1] - 0.3 * shock[i-2] + shock[i]
  Ysim7[i] =   0.6 * shock[i-1] - 0.3 * shock[i-2] + shock[i]
  Ysim8[i] =  -0.6 * shock[i-1] + 0.3 * shock[i-2] + shock[i]
}

par(mfrow=c(2,2))
ts.plot(Ysim5, main = 'positive coeff')
ts.plot(Ysim6, main = 'negative coeff')
ts.plot(Ysim7, main = 'alternate coeff (+, -)')
ts.plot(Ysim8, main = 'alternate coeff(-, +)')


# ************************************************* ******************
##### AUTO-REGRESSIVE PROCESSES #####
# Let's start with AR processes (1)
# ************************************************* ******************
# As with MA processes, we use the 'arima.sim' command
Y1 = arima.sim (n = n, list (ar = 0.9), sd = sqrt (1))
Y1.short = arima.sim (n = n, list (ar = 0.2), sd = sqrt (1)) # short-memory process
# and with negative coefficient
Y2 = arima.sim (n = n, list (ar = -0.9), sd = sqrt (1))

# let's see the effect of the coefficient on the series
par (mfrow = c (2,1))
ts.plot (Y1, main = 'coeff. 0.9')
ts.plot (Y1.short, main = 'coeff. 0.2')

# let's see the effect of the coefficient on the ACF
acf (Y1, main = 'coeff. 0.9')
acf (Y1.short, main = 'coeff. 0.2')

# Compare the acf of the two processes:
# one with a positive coefficient and the other with a negative coefficient
acf (Y1, main = 'positive coefficient')
acf (Y2, main = 'negative coefficient')

# Compare the pacf of the two processes
pacf (Y1, main = 'positive coefficient')
pacf (Y2, main = 'positive coefficient')

# we can get the same processes by simulating them without the command 'arima.sim'
# we initialize the vectors that will contain the data
Ysim1 = c (); Ysim1 [1] = 0
Ysim2 = c (); Ysim2 [1] = 0

for (i in 2: n) {
  Ysim1 [i] = 0.9 * Ysim1 [i-1] + rnorm (1, 0, 1)
  Ysim2 [i] = -0.9 * Ysim2 [i-1] + rnorm (1, 0, 1)
}

par (mfrow = c (2,1))
ts.plot (Ysim1, main = 'positive coefficient')
ts.plot (Ysim2, main = 'negative coefficient')

# Compare the acf of the two processes
acf (Ysim1, main = 'positive coefficient')
acf (Ysim2, main = 'negative coefficient')

# Compare the pacf of the two processes
pacf (Ysim1, main = 'positive coefficient')
pacf (Ysim2, main = 'negative coefficient')

# Stationarity is necessary to study a time series properties.
# Simulate a non-stationary AR (1) process with coefficient 1.1
Y3 = c (); Y3 [1] = 0
for (i in 2: n) {
  Y3 [i] = 1.1 * Y3 [i-1] + rnorm (1, 0, 1)
  
par (mfrow = c (1,1))
ts.plot (Y3)
# the process explodes: it diverges to +/- infinity as n increases

# there is a particular type of non-stationary process called Random Walk, in which
# the coefficient determining the impact of each observation on the next is equal to 1
Y4 = c (); Y4 [1] = 0
for (i in 2: n) {
  Y4 [i] = 1 * Y4 [i-1] + rnorm (1, 0, 1)
}

par (mfrow = c (1,1))
ts.plot (Y4)

# Compare the acf of a stationary process with a Random Walk
par (mfrow = c (2,1))
acf (Y1, main = 'Stationary process')
acf (Y4, main = 'Random Walk process')


# The Dickey-Fuller test (adf.test): the null hypothesis (H0) is that the process is a Random Walk (default)
# NB! The command 'adf.test' allows to choose the alternative hypothesis (H1):
library (Tseries)
adf.test (Y4, alternative = "stationary")
# NB! The D-F test tests only non-stationarity in covariance (unit root)

# Let's move on to AR processes (2), we can have four cases:
# NOTE! See stationarity triangle
# 1. coefficients with alternate signs (-, +)
Y5 = arima.sim (n = n, list (ar = c (-0.6, 0.3)), sd = sqrt (1))
# 2. both positive coefficients
Y6 = arima.sim (n = n, list (ar = c (0.6, 0.3)), sd = sqrt (1))
# 3. both negative coefficients and complex roots
Y7 = arima.sim (n = n, list (ar = c (-1.6, -0.8)), sd = sqrt (1))
# 4. coefficients with alternate signs (+, -) and complex roots
Y8 = arima.sim (n = n, list (ar = c (1.6, -0.8)), sd = sqrt (1))

par (mfrow = c (2,2))
ts.plot (Y5, main = 'alternate signs (-, +)')
ts.plot (Y6, main = 'both positive')
ts.plot (Y7, main = 'negatives and compl. roots')
ts.plot (Y8, main = 'alternate (+, -) and compl. roots')

# Compare the acf of the four processes
acf (Y5, main = 'alternate signs (-, +)')
acf (Y6, main = 'both positive')
acf (Y7, main = 'negatives and compl. roots')
acf (Y8, main = 'alternate (+, -) and compl. roots')

# Compare the pacf of the four processes
pacf (Y5, main = 'alternate signs (-, +)')
pacf (Y6, main = 'both positive')
pacf (Y7, main = 'negatives and compl. roots')
pacf (Y8, main = 'alternate (+, -) and compl. roots')

# the same processes are obtained by simulating them without the command 'arima.sim':
Ysim5 = c(); Ysim5[1] = 0; Ysim5[2] = 0; 
Ysim6 = c(); Ysim6[1] = 0; Ysim6[2] = 0; 
Ysim7 = c(); Ysim7[1] = 0; Ysim7[2] = 0; 
Ysim8 = c(); Ysim8[1] = 0; Ysim8[2] = 0; 

for (i in 3:n){
  Ysim5[i] =  -0.6 * Ysim5[i-1] + 0.3 * Ysim5[i-2] + rnorm(1, 0, 1)
  Ysim6[i] =   0.6 * Ysim6[i-1] + 0.3 * Ysim6[i-2] + rnorm(1, 0, 1)
  Ysim7[i] =  -1.6 * Ysim7[i-1] - 0.8 * Ysim7[i-2] + rnorm(1, 0, 1)
  Ysim8[i] =   1.6 * Ysim8[i-1] - 0.8 * Ysim8[i-2] + rnorm(1, 0, 1)
}
par (mfrow = c (2,2))
ts.plot (Ysim5, main = 'alternate signs (-, +)')
ts.plot (Ysim6, main = 'both positive')
ts.plot (Ysim7, main = 'negatives and compl. roots')
ts.plot (Ysim8, main = 'alternate (+, -) and compl. roots')


##### EXAMPLES OF AR PROCESSES #####
library (Tseries)
date (NelPlo)
? NelPlo

# Consider the time series of the consumer price index from 1860 to 1988 (annual frequency)
par (mfrow = c (3,1))
plot (cpi)
acf (cpi) # makes us think of a very persistent autoregressive process
PACF (cpi)

# How could the similar simulated process be? Try it yourself !!!
n = length (cpi)
YsimSpec = c ()
# ...
# for (i in ??: n) {
# YsimSpec [i] = ...
#}

ts.plot (YsimSpec)
acf (YsimSpec)
PACF (YsimSpec)

# Another example
plot (unemp)
acf (unemp) # Could be an AR (2) with complex roots
pacf (unemp) # with alternate signs (+, -) such as the case Y8