(a)
# Obtain the monthly time series
library(tidyquant)
## Loading required package: lubridate
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
##
## date
## Loading required package: PerformanceAnalytics
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
## Loading required package: quantmod
## Loading required package: TTR
## Version 0.4-0 included new data defaults. See ?getSymbols.
## Loading required package: tidyverse
## Warning: package 'tidyverse' was built under R version 3.4.4
## -- Attaching packages ----------------------------------------------------------- tidyverse 1.2.1 --
## v ggplot2 2.2.1 v purrr 0.2.4
## v tibble 1.4.2 v dplyr 0.7.4
## v tidyr 0.8.0 v stringr 1.3.0
## v readr 1.1.1 v forcats 0.3.0
## Warning: package 'ggplot2' was built under R version 3.4.4
## -- Conflicts -------------------------------------------------------------- tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks xts::first()
## x lubridate::intersect() masks base::intersect()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks xts::last()
## x lubridate::setdiff() masks base::setdiff()
## x lubridate::union() masks base::union()
library(KFAS)
## Warning: package 'KFAS' was built under R version 3.4.4
data1 <-
tq_get("TB3MS", get = "economic.data",
from = "1980-01-01", to = "2017-12-31") %>%
rename(TB3MS = price) %>%
mutate(lTB3MS = log(TB3MS))
data2 <-
tq_get("^GSPC", get = "stock.prices", collapse = "monthly",
from = "1980-01-01", to = "2017-12-31") %>%
select(date, adjusted) %>%
mutate(monyear = as.yearmon(date)) %>%
group_by(monyear) %>%
summarise(SP500 = mean(adjusted)) %>%
ungroup()
data3 <-
tq_get("IBM", get = "stock.prices", collapse = "monthly",
from = "1980-01-01", to = "2017-12-31") %>%
select(date, adjusted) %>%
mutate(monyear = as.yearmon(date)) %>%
group_by(monyear) %>%
summarise(IBM = mean(adjusted)) %>%
ungroup()
(b)
# Construct the time series
rSP500 <- ts(100*diff(log(data2$SP500), differences = 1), start=c(1980,1), frequency = 12)
rIBM <- ts(100*diff(log(data3$IBM), differences = 1), start=c(1980,1), frequency = 12)
rTB3MS <- ts((data1$TB3MS)/12, start=c(1980,1), frequency = 12)
erSP500 <- rSP500 - rTB3MS # Excess Return for S&P500
erIBM <- rIBM - rTB3MS # Excess Return for IBM
par(mfrow=c(3,1))
plot(rSP500, main=expression("Monthly S&P500 Returns"))
plot(rIBM, main=expression("Monthly IBM Returns"))
plot(rTB3MS, main=expression("3 Month Tbill"))
par(mfrow=c(2,1))
plot(erSP500, main=expression("Excess Return for S&P500"))
plot(erIBM, main=expression("Excess Return for IBM"))
(c)
# Estimating CAPM - OLS
# Estimate a relating excess returns for IBM and excess returns for S&P500
# Since series do not align in time, we restrict them
# over the window from Feb 1980 to Dec 2017.
rSP5002 <- window(rSP500, start = 1980 + 2/12, end = 2017 + 12/12)
## Warning in window.default(x, ...): 'end' value not changed
rIBM2 <- window(rIBM, start = 1980 + 2/12, end = 2017 + 12/12)
## Warning in window.default(x, ...): 'end' value not changed
rTB3MS2 <- window(rTB3MS, start = 1980 + 2/12, end = 2017 + 12/12)
## Warning in window.default(x, ...): 'end' value not changed
erSP5002 <- rSP5002 - rTB3MS2 # Excess Return for SP500
erIBM2 <- rIBM2 - rTB3MS2 # Excess Return for IBM
OLS <- lm(erIBM2 ~ erSP5002)
alpha.ols <- OLS$coefficients[1]
beta.ols <- OLS$coefficients[2]
plot(erIBM2,erSP5002, main="OLS")
abline(lm(erSP5002~erIBM2))
summary(OLS)
##
## Call:
## lm(formula = erIBM2 ~ erSP5002)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.5289 -2.6650 0.1015 2.6284 19.0401
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.28485 0.23545 1.21 0.227
## erSP5002 0.94846 0.06563 14.45 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.987 on 451 degrees of freedom
## Multiple R-squared: 0.3165, Adjusted R-squared: 0.315
## F-statistic: 208.9 on 1 and 451 DF, p-value: < 2.2e-16
(d)
# Estimating CAPM - KFS
# Estimate a CAPM model with time varying coefficients at and bt using Kalman filter.
# Create a plot showing how the smoothed states at and bt changed over time.
# SS stands for State-Space
SS.erSP5002 <- as.ts(erSP5002)
SS.erIBM2 <- as.ts(erIBM2)
# get number of observatons
T <- length(SS.erSP5002)
# construct system matrices for state-space model - CAPM with time variables alpha and beta
Zt <- array(rbind(rep(1,T), SS.erSP5002), dim=c(1,2,T))
Ht <- matrix(NA)
Tt <- diag(2)
Rt <- diag(2)
Qt <- matrix(c(NA,0,0,NA), 2,2)
# use diffuse prior
a1 <- matrix(c(0, 1), 2, 1)
P1 <- matrix(0, 2, 2)
P1inf <- diag(2)
# define state-space CAPM model
y.SS <- SSModel(SS.erIBM2 ~ -1 + SSMcustom(Z=Zt, T=Tt, R=Rt, Q=Qt, a1=a1, P1=P1, P1inf=P1inf), H=Ht)
# estimate variances of innovations using maximum likelihood
y.SS.ML <- fitSSM(y.SS, inits = c(0.001,0.001,0.001), method="BFGS")
# Kalman filtering and smoothing, with parameters in Q and H set to maximum likelihood estimates
y.SS.KFS <- KFS(y.SS.ML$model)
# To extract the time-varying coefficients from this model we use
alpha.KFS <- ts( cbind(y.SS.KFS$a[,1], y.SS.KFS$alphahat[,1]), start=c(1980,2), frequency=12)
beta.KFS <- ts( cbind(y.SS.KFS$a[,2], y.SS.KFS$alphahat[,2]), start=c(1980,2), frequency=12)
par(mfrow=c(2,1), mar=c(2,2,2,1))
plot.ts(alpha.KFS, plot.type="single", xlab="",ylab="Alpha", col=c("blue","red"), lwd=2)
legend("topright", c("filtered","smoothed"), col=c("blue","red"), lwd=2, cex=0.7, bty="n")
plot.ts(beta.KFS, plot.type="single", xlab="",ylab="Beta", col=c("blue","red"), lwd=2)
legend("topright", c("filtered","smoothed"), col=c("blue","red"), lwd=2, cex=0.7, bty="n")
par(mfrow=c(2,1), mar=c(2,2,2,1))
plot.ts(alpha.KFS[,2], plot.type="single", xlab="",ylab="Alpha", col="red", lwd=2)
legend("topright", "smoothed", col="red", lwd=2, cex=0.7, bty="n")
plot.ts(beta.KFS[,2], plot.type="single", xlab="",ylab="Beta", col="red", lwd=2)
legend("topright", "smoothed", col="red", lwd=2, cex=0.7, bty="n")
(e)
# Compare at and bt
mean(OLS$coefficients[1])
## [1] 0.2848487
mean(OLS$coefficients[2])
## [1] 0.9484601
mean(y.SS.KFS$a[, 1])
## [1] 0.2869122
mean(y.SS.KFS$a[, 2])
## [1] 0.9919725
mean(y.SS.KFS$alphahat[, 1])
## [1] 0.2958768
mean(y.SS.KFS$alphahat[, 2])
## [1] 0.9492199
# The coefficients of KFS and OLS estimates are nearly identical; however, the alpha and beta from the OLS estimate is smaller than the KFS estimate.