1 General directions for this Workshop

You will work in RStudio. Create an R Notebook document (File -> New File -> R Notebook), where you have to write whatever is asked in this workshop.

At the beginning of the R Notebook write Workshop 2 - Financial Econometrics II and your name (as we did in previous workshop).

Create an R Notebook document (File -> New File -> R Notebook), where you have to write whatever is asked in this workshop. More specifically, you have to:

  • Replicate all the R Code along with its output.

  • You have to do whatever is asked in the workshop. It can be: a) Responses to specific questions and/or do an exercise/challenge.

Any QUESTION or any STEP you need to do will be written in CAPITAL LETTERS. For ANY QUESTION, you have to RESPOND IN CAPITAL LETTERS right after the question.

  • It is STRONGLY RECOMMENDED that you write your OWN NOTES as if this were your personal notebook. Your own workshop/notebook will be very helpful for your further study.

You have to keep saving your .Rmd file, and ONLY SUBMIT the .html version of your .Rmd file. Pay attention in class to know how to generate an html file from your .Rmd.

2 Non stationary variables - The Random Walk model for stock prices

The random walk hypothesis in Finance (Fama, 1965) states that the natural logarithm of stock prices behaves like a random walk with a drift. A random walk is a series (or variable) that cannot be predicted. Imagine that \(Y_t\) is the log price of a stock for today (t). The value of Y for tomorrow (\(Y_{t+1}\)) will be equal to its today’s value (\(Y_t\)) plus a constant value (\(φ_0\)) plus a random shock. This shock is a pure random value that follows a normal distribution with mean=0 and a specific standard deviation \(σ_ε\). The process is supposed to be the same for all future periods. In mathematical terms, the random walk model is the following:

\[ Y_t = φ_0 + Y_{t−1} + ε_t \]

The \(ε_t\) is a random shock for each day, which is the result of the log price movement due to all news (external and internal to the stock) that influence the price. \(φ_0\) refers as the drift of the series. If \(|φ_0|\) > 0 we say that the series is a random walk with a drift. If \(φ_0\) is positive, then the variable will have a positive trend over time; if it is negative, the series will have a negative trend.

If we want to simulate a random walk, we need the values of the following parameters/variables:

  • \(Y_0\), the first value of the series
  • \(φ_0\), the drift of the series
  • \(σ_ε\), the standard deviation (volatility) of the random shock

2.1 Q Monte Carlo simulation for the random walk model

Let’s go and run a Monte Carlo simulation for a random walk of the S&P 500. We will use real values of the S&P500 to estimate the previous 3 parameters.

2.1.1 Loading/installing R packages

The following R packages need to be installed for most of the class workshops:

  • fpp2

  • fpp3

  • quantmod

  • dplyr

  • ggplot2

Go to the right-bottom windows of RStudio, select the Package tab, click install, and install both R packages. These packages include many other R packages for time-series data management and analysis.

These fpp2 and fpp3 packages were written by Rob J Hyndman and George Athanasopoulos, professors from Monash University at Australia. They are also business consultants with many years of experience doing both serious research in time-series and also applying their findings in the real world.

The quantmod package was written by Jeffrey A. Ryan, one of the most important contributors to R packages for Finance, and organizer of the famous R/Finance Conference held in Chicago each year since 2009.

Once you install these packages, load them in memory:

library(quantmod)
library(fpp3)
library(dplyr)
library(ggplot2)
library(zoo)

2.1.2 Downloading data for the S&P500

Download the S&P500 historical daily data from Yahoo Finance from 2009 to date.

getSymbols("^GSPC", from="2009-01-01")
## [1] "^GSPC"

Now we generate the log of the S&P index using the closing price/quotation, and create a variable N for the number of days in the dataset:

lnsp<-log(Ad(GSPC))
# I assign a name for the index:
names(lnsp)<-c("lnsp")
N<-nrow(lnsp)

Now we will simulate 2 random walk series estimating the 3 parameters from this log series of the S&P500:

  1. a random walk with a drift (name it rw1), and

  2. a random walk with no drift (name it rw2).

2.1.3 Estimating the parameters of the random walk model

We have to consider the mathematical definition of a random walk and estimate its parameters (initial value, phi0, volatility of the random shock) from the real daily S&P500 data.

Now, we create a variable for a random walk with a drift trying to model the log of the S&P500.

Reviewing the random walk equation again:

\[ Y_t = φ_0 + Y_{t−1} + ε_t \] The \(ε_t\) is the random shock of each day, which represents the overall average perception of all market participants after learning the news of the day (internal and external news announced to the market).

Remember that \(\varepsilon_{t}\) behaves like a random normal distributed variable with mean=0 and with a specific standard deviation \(\sigma_{\varepsilon}\).

For the simulation of the random walk, you need to estimate the values of

  • \(y_{0}\), the first value of the series, which is the log S&P500 index of the first day

  • \(\phi_{0}\)

  • \(\sigma_{\varepsilon}\)

You have to estimate \(\phi_{0}\) using the last and the first real values of the series following the equation of the random walk.

Here you can see possible values of a random walk over time:

\[ Y_{0} = Initial value \] \[ Y_{1} = \phi_{0} + Y_{0} + \varepsilon_{1} \] \[ Y_{2} = \phi_{0} + Y_{1} + \varepsilon_{2} \] Substituting \(Y_{1}\) with its corresponding equation: \[ Y_{2} = \phi_{0} + \phi_{0} + Y_{0} + \varepsilon_{1} + \varepsilon_{2} \] Re-arranging the terms: \[ Y_{2} = 2*\phi_{0} + Y_{0} + \varepsilon_{1} + \varepsilon_{2} \] If you continue doing the same until the last N value, you can get: \[ Y_{N} = N*\phi_{0} + Y_{0} + \sum_{t=1}^{N}\varepsilon_{t} \]

This mathematical result is kind of intuitive. The value of a random walk at time N will be equal to its initial value plus N times phi0 plus the sum of ALL random shocks from 1 to N.

Since the mean of the shocks is assumed to be zero, then the expected value of the sum of the shocks will also be zero. Then:

\[ E[Y_{N}] = N*\phi_{0} + Y_{0} \] From this equation we see that \(phi_{0}\) can be estimated as: \[ \phi_{0} = \frac{(Y_{N} - Y_{0})}{N} \]

Then, \(\phi_{0}\) = (last value - first value) / # of days.

I calculate \(\phi_{0}\) following this formula:

phi0<- (as.numeric(lnsp$lnsp[N])-as.numeric(lnsp$lnsp[1])) / N
cat("The value for phi0 is ",phi0)
## The value for phi0 is  0.0004626235

Remember that N is the total # of observations, so lnsp[N] has last daily value of the log of the S&P500.

Now we need to estimate sigma, which is the standard deviation of the shocks. We can start estimating its variance first. It is known that the variance of a random walk cannot be determined unless we consider a specific number of periods.

Then, let’s consider the equation of the random walk series for the last value (\(Y_N\)), and then estimate its variance from there:

\[ Y_{N} = N*\phi_{0} + Y_{0} + \sum_{t=1}^{N}\varepsilon_{t} \] Using this equation, we calculate the variance of \(Y_N\) :

\[ Var(Y_{N}) = Var(N*\phi_{0}) + Var(Y_{0}) + \sum_{t=1}^{N}Var(\varepsilon_{t}) \] The variance of a constant is zero, so the first two terms are equal to zero.

Now analyze the variance of the shock:

Since it is supposed that the volatility (standard deviation) of the shocks is about the same over time, then:

\[ Var(\varepsilon_{1}) = Var(\varepsilon_{2}) = Var(\varepsilon_{N}) = \sigma_{\varepsilon}^2 \] Then the sum of the variances of all shocks is actually the variance of the shock times N. Then the variance of all the shocks is actually the variance of \(Y_N\).

Then we can write the variance of \(Y_N\) as:

\[ Var(Y_{N}) = N * Var(\varepsilon)= N*\sigma_{\varepsilon}^2 \] To get the standard deviation of \(Y_N\) we take the square root of the variance of \(Y_N\):

\[ SD(Y_{N}) = \sqrt{N}*SD(\varepsilon) \] We use sigma character for standard deviations:

\[ \sigma_{Y} = \sqrt{N}*\sigma_{\varepsilon} \]
Finally we express the volatility of the shock (\(\sigma_{\varepsilon}\)) in terms of the volatility of \(Y_N\) (\(\sigma_{Y}\)):

\[ \sigma_{\varepsilon} = \frac{\sigma_{Y}}{\sqrt{N}} \]

Then we can estimate sigma as: sigma = StDev(lnsp) / sqrt(N). Let’s do it:

sigma<-sd(lnsp$lnsp) / sqrt(N)
cat("The volatility of the log is = ",sd(lnsp$lnsp),"\n")
## The volatility of the log is =  0.4360899
cat("The volatility for the shock is = ",sigma)
## The volatility for the shock is =  0.007582166

2.1.4 Simulating the random walk with drift

Now you are ready to start the simulation of random walk using rw1: \[ rw1_{t} = \phi_{0} + rw1_{t-1} + \varepsilon_{t} \]

The \(\phi_{0}\) coefficient is also drift of the random walk.

We will create a new column in the lnsp R dataset for the random walk with the name rw1.

lnsp$rw1 = 0

I assigned zero to all values before I do the simulation.

I start assigning the first value of the random walk to be equal to the first value of the log of the S&P500:

lnsp$rw1[1]<-lnsp$lnsp[1]

Now assign random values from day 2 to the last day following the random walk. For each day, we create the random shock using the function rnorm. We create this shock with standard deviation equal to the volatility of the shock we calculated above (the sigma). We indicate that the mean =0:

shock <- rnorm(n=N,mean=0,sd=sigma)
lnsp$shock<-shock

We can see the shock over time:

plot(shock, type="l", col="blue")

We can also see whether the shock behaves like a normal distribution by doing its histogram:

hist(lnsp$shock)

As expected, the shock behaves similiar to a normal-distributed variable.

Now we are ready to start the simulation of random walk. Then we fill the values for rw1. Remembering the formula for the random walk process:

\[ rw1_{t} = \phi_{0} + rw1_{t-1} + \varepsilon_{t} \] We start the random walk with the first value of the log of the S&P500. Then, from day 2 we do the simulation according to the previous formula and using the random shock just created:

# I create separate vectors:
rw1<-single(length=N)

# I assign the first value of the random-walk to be equal to real log value of the S&P
rw1[1]<-as.numeric(lnsp$lnsp[1])
# Now from day 1 I generate the values of the random walk following the formula:
for (i in 2:N){
  rw1[i] <- phi0 + rw1[i-1] + shock[i] 
}
lnsp$rw1<-rw1

I plot the simulated random walk and the real log of the S&P500:

ts.plot(lnsp$rw1)
lines(seq(1,N),lnsp$lnsp, col="blue")

2.1.5 Simulating a random walk with no drift

Now we can do a simulation but now without the drift. I this case, the \(\phi_{0}\) coefficient must be zero.

Use another variable rw2 for this. You can follow the logic we did for rw1, but now \(\phi_{0}\) will be equal to zero, so we do not include it into the equation:

rw1_v2<-single(length=N)
rw1_v2[1]<-lnsp$lnsp[1]
for (i in 2:N){
  rw1_v2[i] <- rw1_v2[i-1] + shock[i] 
}

ts.plot(lnsp$lnsp, col="blue")
# I plot both lines to compare 
lines(rw1_v2, col="green")

WHAT DO YOU OBSERVE with this plot? EXPLAIN WITH YOUR WORDS.

Now run a simple regression to check whether the rw1 is statistically related to the log of the S&P500. Use rw1 as explanatory variable. Show the regression results as comments.

regmodel<-lm(lnsp$lnsp~lnsp$rw1)
# I see statistics on my regression model
s_regmodel <- summary(regmodel)
s_regmodel
## 
## Call:
## lm(formula = lnsp$lnsp ~ lnsp$rw1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.57979 -0.05250  0.01284  0.07292  0.21532 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.706303   0.031222   22.62   <2e-16 ***
## lnsp$rw1    0.938743   0.004246  221.07   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1098 on 3306 degrees of freedom
## Multiple R-squared:  0.9366, Adjusted R-squared:  0.9366 
## F-statistic: 4.887e+04 on 1 and 3306 DF,  p-value: < 2.2e-16

DOES THE REGRESSION RESULT MAKE SENSE? EXPLAIN WHY YES OR WHY NOT?

DOES THE LOG OF THE S&P500 LOOKS LIKE A RANDOM WALK? WHY YES OR WHY NOT?

# I plot the natural log pf S&P500
ts.plot(lnsp$lnsp)
# I plot my random walk with a drift
lines(rw1, col="blue")

DO YOU THINK THAT WE CAN USE THIS TYPE OF SIMULATION TO PREDICT STOCK PRICES OR INDEXES? WHY YES OR WHY NOT?

ts.plot(lnsp$lnsp,col="blue")
lines(rw1_v2, col="green")
lines(rw1)

3 SIMULATING A RANDOM WALK AND AN AR(1) PROCESS

3.1 Create the dataset with simulation

  1. Create a dataset of 1000 observations
rm(list = ls())
obs = 1:1000
#View(obs)
  1. Create a Gaussian white noise variable with 1,000 values with variance 0.09. Call this variable e1. The Gaussian white noise is a normal variable with mean=0 and a specic variance. We will use Gaussian white noises to model financial returns.
e1=rnorm(1000,0,sqrt(0.09))
N<-length(obs)
head(e1)
## [1]  0.08558263  0.10620418 -0.39460822  0.26273330 -0.23098207 -0.52738559
tail(e1)
## [1]  0.3694736 -0.3610698  0.1220510  0.4141831  0.2447920 -0.1811947
hist(e1)

3.2 Simulating an AR(1) with phi0=1 and phi1=0.7

  1. Following the formula for a Simple or First-order Autoregressive, AR(1):

  1. Use the Gaussian white noise created (e1) for this

  2. Declare a variable for the coefficient φ1(phi1) and equal this coefficient to 0.7.

  3. Declare a variable for the coefficient φ0(phi0) and equal this coefficient to 1

  4. Using simple simulation, generate the variable y1 as an AR(1) model using the above terms

  5. The variable y1 is an AR(1) process or model. Graph this variable over time

You can run the following code for the previous steps:

e1=rnorm(1000,0,sqrt(0.09))

y1<-single(length=N)
phi1 <- 0.7
phi0 <- 1

# Now from day 1 I generate the values of the random walk following the formula:
for (i in 2:N){
  y1[i] <- phi0 + phi1*y1[i-1] + e1[i-1] 
}

ts.plot(y1, col="darkblue")

a). WHAT DO YOU SEE? LOOKING AT THE GRAPH, DOES THE MEAN OF THE SERIES CONVERGE TO A VALUE? IF YES, WHICH VALUE?

b). WHAT IS THE EXPECTED VALUE OF y1 ACCORDING TO THE AR(1) MODEL? PROVIDE THE FORMULA AND CALCULATE THE EXPECTED VALUE.

c). IS THE EXPECTED VALUE OF y1 SIMILAR TO THE MEAN YOU SAW IN THE FIRST GRAPH?

d). IS THE VOLATILITY (STANDARD DEVIATION) SIMILAR IN ALL TIME PERIODS?

3.3 Simulating an AR(1) with phi0=0 and phi1=0.7

Now generate another series y2 with the same parameters than y1, but just change the constant phi0 from 1 to zero.

  1. Graph y2 over time.

  2. Is this time series also an AR(1)?

  3. WHAT IS THE EXPECTED VALUE OF y2? ACCORDING TO THE AR(1) MODEL? PROVIDE THE FORMULA AND CALCULATE THE EXPECTED VALUE.

e2=rnorm(1000,0,sqrt(0.09))

y2<-single(length=N)
phi1 <- 0.7
phi0 <- 0

# Now from day 1 I generate the values of the random walk following the formula:
for (i in 2:N){
  y2[i] <- phi0 + phi1*y2[i-1] + e2[i-1] 
}

ts.plot(y2, col="darkred")

3.4 Simulating a Random walk with phi0=0

Now generate another series y3, but now modify the parameters to simulate a Random Walk process with phi0=0.

  1. Graph y3 over time. Observe its behaviour.
  2. IS THE MEAN OF THE SERIES CONSTANT OVER TIME?
  3. IS THE VOLATILITY OF THE SERIES CONSTANT OVER TIME? EXPLAIN
e3=rnorm(1000,0,sqrt(0.09))
y3<-single(length=N)
phi1 <- 1
phi0 <- 0

# Now from day 1 I generate the values of the random walk following the formula:
for (i in 2:N){
  y3[i] <- phi0 + phi1*y3[i-1] + e3[i-1] 
}

ts.plot(y3, col="yellow")

3.5 Simulating an AR(1) with phi0=0 and phi1=0.99

Now generate another series y4 as an AR(1) process with phi0=0 and phi1=0.99

  1. Graph y4 over time.

  2. IS THIS SERIES AN AR(1) ? EXPLAIN WHY YES OR WHY NOT

  3. DOES THE MEAN CONVERGE TO A SPECIFIC VALUE? IF YES, TO WHICH ONE?

e4=rnorm(1000,0,sqrt(0.09))

y4<-single(length=N)
phi1 <- 0.99
phi0 <- 0

# Now from day 1 I generate the values of the random walk following the formula:
for (i in 2:N){
  y4[i] <- phi0 + phi1*y4[i-1] + e4[i-1] 
}

ts.plot(y4, col="orange")

3.6 Weakly stationarity.

A time series is weakly stationary if:

  • its expected value is constant over time (about the same over time)

  • its expected variance over time is constant

  • the covariance (or correlation) between y and y(t+h) is the same for any t and any h

3.6.1 Checking for weakly stationary:

We will use the rollaply function from the zoo package to calculate rolling windows and estimate rolling means and rolling standard deviations.

  1. Using the rollaply function compute the rolling means using 12 periods for the series y1, y2, y3, and y4:
y1_mean<-rollapply(y1,12,mean)
y2_mean<-rollapply(y2,12,mean)
y3_mean<-rollapply(y3,12,mean)
y4_mean<-rollapply(y4,12,mean)

These new datasets will have 989 periods instead of 1,000 since we calculated rolling means so that we calculate the first mean until the row 12 for each variable.

  1. Merge the rolling means into one data frame to better manage the time series:
# Create a variable for the time period:
x <- seq(1,989,1)
# Merge the x and the rolling means datasets into one dataset:
all_means <- tbl_df(data.frame(x, y1_mean, y2_mean, y3_mean, y4_mean))
## Warning: `tbl_df()` was deprecated in dplyr 1.0.0.
## Please use `tibble::as_tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
  1. Graph the 12-period rolling means to see how the mean moves over time:
p <- ggplot(all_means, aes(x = x))

p+ geom_line(aes(y=y1_mean), colour="darkblue")+
  geom_line(aes(y=y2_mean), color="darkred")+
  geom_line(aes(y=y3_mean), color="yellow")+
  geom_line(aes(y=y4_mean), color="orange")

  1. Using the rollaply function and a window of 12 time periods, compute the rolling standard deviations of the series y1, y2, y3, and y4:
y1_sd<-rollapply(y1,12,sd)
y2_sd<-rollapply(y2,12,sd)
y3_sd<-rollapply(y3,12,sd)
y4_sd<-rollapply(y4,12,sd)
  1. Graph the rolling standard deviations of each series to see how the standard deviation moves for each window:
all_sd <- tbl_df(data.frame(x, y1_sd, y2_sd, y3_sd, y4_sd))
pl <- ggplot(data = all_sd, aes(x = x))

pl+ geom_line(aes(y=y1_sd), color="darkblue")+
  geom_line(aes(y=y2_sd), color="darkred")+
  geom_line(aes(y=y4_sd), color="yellow")+
  geom_line(aes(y=y3_sd), color="orange")

  1. WHICH OF THE SERIES (y1, y2, y3, y4) IS (ARE) WEAKLYSTATIONARY AND WHICH IS (ARE) NOT? BRIEFLY EXPLAIN

4 Reading

Go to the course site and carefully re-read the Note Basics of time-series au- toregressive models and read the Note: Introduction to ARMA Models.

5 Quiz 2 and W2 submission

Go to Canvas and respond Quiz 2. You will have 3 attempts. Questions in this Quiz are related to concepts of the readings related to this Workshop.

The grade of this Workshop will be the following:

  • Complete (100%): If you submit an ORIGINAL and COMPLETE HTML file with all the activities, with your notes, and with your OWN RESPONSES to questions
  • Incomplete (75%): If you submit an ORIGINAL HTML file with ALL the activities but you did NOT RESPOND to the questions and/or you did not do all activities and respond to some of the questions.
  • Very Incomplete (10%-70%): If you complete from 10% to 75% of the workshop or you completed more but parts of your work is a copy-paste from other workshops.
  • Not submitted (0%)

Remember that you have to submit your .html file through Canvas BEFORE NEXT CLASS.