Workshop 5
Andrea Contel
Download DAILY stock data of Gruma (GRUMAB.MX) and the IPyC (^MXX) from Yahoo from 2017 to March 11, 2020:
# Remember to load the library you will be using
library(quantmod)## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo# Obtain daily data
getSymbols(c("GRUMAB.MX", "^MXX"), from="2017-01-01",
to="2020-03-11", periodicity="daily", src="yahoo")## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.## [1] "GRUMAB.MX" "^MXX"# Visualize data
head(GRUMAB.MX, 5)## GRUMAB.MX.Open GRUMAB.MX.High GRUMAB.MX.Low GRUMAB.MX.Close
## 2017-01-02 263.50 264.00 262.01 263.10
## 2017-01-03 263.97 265.00 257.50 258.77
## 2017-01-04 260.06 266.06 257.00 265.74
## 2017-01-05 265.74 275.00 264.46 274.41
## 2017-01-06 266.00 275.00 266.00 274.64
## GRUMAB.MX.Volume GRUMAB.MX.Adjusted
## 2017-01-02 22518 240.7950
## 2017-01-03 420265 236.8321
## 2017-01-04 1120534 243.2112
## 2017-01-05 1175996 251.6196
## 2017-01-06 860582 251.8306head(MXX, 5)## MXX.Open MXX.High MXX.Low MXX.Close MXX.Volume MXX.Adjusted
## 2017-01-02 45642.80 45792.46 45593.26 45695.10 6765800 45695.10
## 2017-01-03 45698.93 46444.72 45690.46 46123.36 164938800 46123.36
## 2017-01-04 46123.36 46587.74 46068.83 46587.74 215385900 46587.74
## 2017-01-05 46586.03 46767.37 46285.88 46719.99 237288300 46719.99
## 2017-01-06 46721.70 46731.90 46027.10 46071.57 201371200 46071.57Visualize the daily stock prices for Gruma:
plot(GRUMAB.MX$GRUMAB.MX.Adjusted)
Now visualize the IPCyC market index:
plot(MXX$MXX.Adjusted)
We can see that Gruma price and the IPCyC index have a negative trend over time, but looking at March 9 it seems that both had important decline.
We can visualize daily returns of Gruma:
r_GRUMA <- diff(log(Ad(GRUMAB.MX)))
plot(r_GRUMA)
We can visualize daily returns of the Market:
r_MXX <- diff(log(Ad(MXX)))
plot(r_MXX)
We can see that the Mexican Market had its worst day at the end of the series, that is March 9th, 2020. Gruma had also a negative return on the same day, but not as dramatic as in the case of the Mexican market index.
Now we can visualize the linear relationship between Gruma returns and the Mexican index return:
plot(as.numeric(r_MXX), as.numeric(r_GRUMA), type = "p",
xlab = "r_MXX", ylab = "r_GRUMA",
main = "Scatter plot between Gruma returns and Market returns")
abline(lm(r_GRUMA~r_MXX), col="red")
Although the scales of both axis have different range values, we can appreciate that the slope of the regression line is not quite big, and it should be less than 1 since the inclination of the line is less than 45 degrees.
Now you have to run the Market Regression Model for Gruma restricting dates up to March 6, 2020 (before the BLACK MONDAY):
# Subset orginal time series of GRUMA returns and assign it to r_GRUMA_bef
r_GRUMA_bef <- r_GRUMA["2017-01-02/2020-03-06"]
# Another way to subset that will give the exact same result
ie <- r_GRUMA[1:799,]
# Do the same for the returns of the market
r_MXX_bef <- r_MXX["2017-01-02/2020-03-06"]
# Use the lm() function to make a linear regression
# bef_bm stands for before black monday
bef_bm <- lm(r_GRUMA_bef ~ as.vector(r_MXX_bef))
# Notice the class of the new object
class(bef_bm)## [1] "lm"# Take a look at the contents of the object
names(bef_bm)## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "na.action" "xlevels" "call" "terms"
## [13] "model"# Use the function summary() to access the names of bef_bm
summ_bef_bm <- summary(bef_bm)
# We are interested in the coefficients of the model
summ_bef_bm$coefficients## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0002884309 0.0004395498 -0.6561962 5.118876e-01
## as.vector(r_MXX_bef) 0.5949982876 0.0500974110 11.8768271 4.567387e-30What you can say about the market risk of GRUMA? EXPLAIN.
#THE MARKET RISK OF GRUMA IS LOW, WE CAN SAY THIS BECAUSE THE STANDARD ERRORS ARE SMALL.
##What is the regression equation of this model? WRITE it as comments in your R script.
gruma_ret <- diff(log(GRUMAB.MX$GRUMAB.MX.Adjusted))
market_ret <- diff(log(MXX$MXX.Adjusted))
reg2 <- lm(gruma_ret ~ market_ret)
summary(reg2)##
## Call:
## lm(formula = gruma_ret ~ market_ret)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043136 -0.007148 -0.000481 0.007422 0.067089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0002807 0.0004386 -0.64 0.522
## market_ret 0.5848178 0.0481871 12.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0124 on 798 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1558, Adjusted R-squared: 0.1548
## F-statistic: 147.3 on 1 and 798 DF, p-value: < 2.2e-16reg2 <- lm(GRUMAB.MX ~ MXX)E(gruma_ret)= b0 + b1 = gruma_ret
E(gruma_ret)= -0.0002 + 0.548 = gruma_ret
Note that your regression DID NOT include the BLACK MONDAY March 9, 2020. Now check what was the real market daily return on March 9th and March 10th, 2020. You can just show the last 5 days of the data set to see the last daily returns:
tail(r_MXX, 5)## MXX.Adjusted
## 2020-03-04 0.02171823
## 2020-03-05 -0.02472212
## 2020-03-06 -0.02283717
## 2020-03-09 -0.06638094
## 2020-03-10 0.02132713tail(r_GRUMA, 5)## GRUMAB.MX.Adjusted
## 2020-03-04 0.020391944
## 2020-03-05 -0.001206561
## 2020-03-06 -0.025035198
## 2020-03-09 -0.030216820
## 2020-03-10 0.010485285Now using the real data of the market return for the BLACK MONDAY, ESTIMATE what would be the expected return for GRUMA on Monday March 9th, 2020 using the regression equation. How GOOD was your prediction compared with the real Gruma return on that day? DO THE CALCULATION of your PREDICTION in your r file.
gruma_ret <- diff(log(GRUMAB.MX$GRUMAB.MX.Adjusted))
market_ret <- diff(log(MXX$MXX.Adjusted))reg2 <- lm(gruma_ret ~ market_ret)
summary(reg2)##
## Call:
## lm(formula = gruma_ret ~ market_ret)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043136 -0.007148 -0.000481 0.007422 0.067089
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0002807 0.0004386 -0.64 0.522
## market_ret 0.5848178 0.0481871 12.14 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0124 on 798 degrees of freedom
## (1 observation deleted due to missingness)
## Multiple R-squared: 0.1558, Adjusted R-squared: 0.1548
## F-statistic: 147.3 on 1 and 798 DF, p-value: < 2.2e-16reg2 <- lm(GRUMAB.MX ~ MXX)Now, use the predict() function to do the prediction of the GRUMA daily returns. Make sure that your manual estimation of the prediction for the BLACK MONDAY is exactly the same you got using the predict() function. To run the predict() function, just type:
predict(bef_bm, newdata = data.frame(r_MXX_bef = -0.06638094))## 1
## -0.03978498You can also create a data frame that contains the values of the explanatory variable you want to try. In order not to copy and paste the value, you can refer to it by subsetting the original time series of the explanatory variable (in this case, r_MXX) However, the explanatory variable in both lm() and predict() must be the same. Thus, in this case, we must convert the subset of r_MXX (an xts object) to a simple numeric value by using as.numeric().
explanatory_var <- data.frame(r_MXX_bef = as.numeric(r_MXX["2020-03-09"]))
# Then, we can use predict() or predict.lm()
# with the explanatory variable object and get the same result.
# ALSO CALLED INDEPENDENT VARIABLE.
predict.lm(bef_bm, newdata = explanatory_var)## 1
## -0.03978498According to the model, GRUMA would have had -3.97% of returns on March 9, 2020 given the real market return of March 9,2020 (-6.63%)
Looking at the real return of GRUMA we displayed above (-3.02%), we can see that the prediction was close to the actual GRUMA return on the Black Monday March 9, 2020 (-3.97% vs -3.02%).
#4 Estimating “moving” betas of the Market Model Systematic risk of a stock is mainly due to movements in the market. Non-systematic risk of a stock is the idiosyncratic risk of the stock that is NOT related to the market. Beta1 represents the market risk of a stock, or the systematic risk of the stock.
We have to be aware that the beta1 and beta0 regression coefficients of a market model of a stock CAN MOVE over time. This is the reason why we always see or calculate the standard error of the beta0 and beta1 in a regression. Remember that the standard error in the context of a regression is the standard deviation of the corresponding beta.
Then, how the estimated betas of a stock move over time? Is the beta1 of a stock stable or does it radically change over time?
We can use the roll_lm() function of the roll package to automatically run many regression models, one for different rolling time window. For example, to run many regressions with a 24-month rolling window for the Gruma monthly returns starting in 2010, you can do the following in R:
We get the monthly data of Gruma from 2010 to March, 2020:
# Start with a clean environment
rm(list=ls())
# Use getSymbols() to download data
getSymbols(c("GRUMAB.MX", "^MXX"), from="2010-01-01",
to="2020-03-31", periodicity="monthly", src="yahoo")## [1] "GRUMAB.MX" "^MXX"# Calculate returns
r_GRUMA <- diff(log(Ad(GRUMAB.MX)))
r_MXX <- diff(log(Ad(MXX)))Now we run the market regression model using windows of 24-month starting in January 2010. Then, R will run many regression models; one regression model for each 24-month window. If I have 10 years and 3 months of data, then I have about 123 months in my data set, so R will run 100 (123-24+1) regression models since there are 100 24-month windows from Jan 2010 to March 2020.
We need to install the roll R package (install it from RStudio/Packages). We will use the roll_lm() function of thie package, indicating the dependent and independent variables of the regression model (in that order: x, y) as vectors or matrices, as well as size of the windows. We can store the resulting data in a new object:
# Load library
library(roll)
# roll_lm() receives only vectors and matrices, so I coerce the data of returns
rollmktmdl <- roll_lm(as.vector(r_MXX), as.vector(r_GRUMA), width = 24)
# Look at the class of the resulting object
class(rollmktmdl)## [1] "list"# Notice the elements of the list object
names(rollmktmdl)## [1] "coefficients" "r.squared" "std.error"As we are interested in the coefficients, we will create an object that contains them:
rolling_betas <- rollmktmdl$coefficients
# The columns have weird names. Let's change them
colnames(rolling_betas) <- c("b0", "b1")If you use View() to take a look at the coefficients, you will notice that the first 24 rows contain NA values. This is because the first 24 rows were taken in the first window and there was no prior information to compute the regression model for these values.
The roll_lm() function returns the coefficients as vectors. However, we are working with time series, so we have to change the class of the coefficients object (rolling_betas).
rolling_betas <- xts(rolling_betas, order.by = index(r_GRUMA))
# Now rolling_betas is an xts objectI graph the beta1 of Gruma for the 100 24-month windows (windows ending in 2012 to 2020):
plot(rolling_betas$b1, main = "Moving beta1 for Gruma stock", ylab = "Beta 1")
##You have to INTERPRET THIS GRAPH (respond in your R script) # FROM THIS GRAPH WE CAN TELL THAT THE BETA1 HAD ITS PEAK FROM 2012 TO 2013 AND AGAIN IN MID 2018, AND THERE WERE GREAT DECLINES. I can also see how beta0 moves over time in the 100 24-month windows:
plot(rolling_betas$b0, main = "Moving beta0 for Gruma stock", ylab = "Beta 0")
##You have to INTERPRET THIS GRAPH (respond in your R script) # FROM THIS GRAPH WE CAN TELL THAT THE GREATEST RETURNS FOT BETA0 ARE FROM 2014 MID 2015, IT FOLLOWS A NORMAL DISTRIBUTION.