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
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(plyr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following objects are masked from 'package:xts':
##
## first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(PerformanceAnalytics)
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
library(tidyverse)
## -- Attaching packages -------------- tidyverse 1.3.0 --
## v ggplot2 3.3.0 v purrr 0.3.3
## v tibble 2.1.3 v stringr 1.4.0
## v tidyr 1.0.0 v forcats 0.4.0
## v readr 1.3.1
## -- Conflicts ----------------- tidyverse_conflicts() --
## x dplyr::arrange() masks plyr::arrange()
## x purrr::compact() masks plyr::compact()
## x dplyr::count() masks plyr::count()
## x dplyr::failwith() masks plyr::failwith()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks xts::first()
## x dplyr::id() masks plyr::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks xts::last()
## x dplyr::mutate() masks plyr::mutate()
## x dplyr::rename() masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(tidyr)
tickers <- c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD")
n <- length(tickers)
data <- new.env()
getSymbols(tickers, src="yahoo", from = "2010-01-01", to = "2018-12-30", env = data, auto.assign = T)
## '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.
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## [1] "SPY" "QQQ" "EEM" "IWM" "EFA" "TLT" "IYR" "GLD"
con = gzcon(url('http://www.systematicportfolio.com/sit.gz','rb'))
source(con)
close(con)
bt.prep(data, align = 'remove.na', fill.gaps = T)
names(data)
## [1] "GLD" "prices" "IWM"
## [4] "dates" "IYR" "TLT"
## [7] "weight" "QQQ" ".getSymbols"
## [10] "symbolnames" "execution.price" "EEM"
## [13] "SPY" "EFA"
prices <- data$prices
n <- ncol(prices)
ret <- prices/mlag(prices) - 1
ret.sample = ret[-1, ]
head(ret.sample)
## EEM EFA GLD IWM
## 2010-01-05 0.007258277 0.0008813503 -0.0009108014 -0.0034386058
## 2010-01-06 0.002092073 0.0042268581 0.0164995902 -0.0009409818
## 2010-01-07 -0.005799118 -0.0038583129 -0.0061878037 0.0073782887
## 2010-01-08 0.007932804 0.0079225530 0.0049630301 0.0054542467
## 2010-01-11 -0.002083333 0.0082096245 0.0132889913 -0.0040296809
## 2010-01-12 -0.016005636 -0.0117810116 -0.0209127164 -0.0108932614
## IYR QQQ SPY TLT
## 2010-01-05 0.002401157 0.0000000000 0.002647093 0.0064580894
## 2010-01-06 -0.000435453 -0.0060318615 0.000704057 -0.0133864256
## 2010-01-07 0.008932440 0.0006501734 0.004221291 0.0016820139
## 2010-01-08 -0.006694040 0.0082304747 0.003327769 -0.0004477891
## 2010-01-11 0.004782630 -0.0040815898 0.001396552 -0.0054877141
## 2010-01-12 -0.016875854 -0.0125108280 -0.009326235 0.0171170828
Q1
mean(ret.sample)
## [1] 0.0002573027
median(ret.sample)
## [1] 0.0006250431
quantile(ret.sample)
## 0% 25% 50% 75% 100%
## -0.0878082616 -0.0050807725 0.0006250431 0.0061394001 0.0925586359
sd(ret.sample)
## [1] 0.0110936
skewness(ret.sample)
## EEM EFA GLD IWM IYR QQQ
## Skewness -0.1802524 -0.4722765 -0.4821824 -0.2216965 -0.1176579 -0.295904
## SPY TLT
## Skewness -0.3951384 -0.1026327
kurtosis(ret.sample)
## EEM EFA GLD IWM IYR QQQ
## Excess Kurtosis 2.789026 5.568187 5.296673 3.740159 6.318948 3.109242
## SPY TLT
## Excess Kurtosis 4.312212 1.556317
Q2
avg.ret <- colMeans(coredata(ret.sample))
cov.d <- cov(coredata(ret.sample))
ones.vec <- matrix(rep(1,n), ncol = 1)
numerator <- solve(cov.d) %*% ones.vec
denominator <- t(ones.vec)%*%numerator
mvp.d <- numerator / as.numeric(denominator)
mvp.d
## [,1]
## EEM -0.12060731
## EFA 0.05284383
## GLD 0.13170178
## IWM -0.05588057
## IYR -0.09957534
## QQQ -0.19637470
## SPY 0.84908458
## TLT 0.43880774
Q3
ret.sample1 <- ret['2010-01-01/2013-12-30', ]
ret.sample.w <- ret.sample1 %>% to.weekly(indexAt = "lastof", OHLC = FALSE)
## Warning in to.period(x, "weeks", name = name, ...): missing values removed
## from data
avg.ret.w <- colMeans(coredata(ret.sample.w))
cov.d.w <- cov(coredata(ret.sample.w))
numerator.w <- solve(cov.d.w) %*% ones.vec
denominator.w <- t(ones.vec)%*%numerator.w
mvp.w <- numerator.w / as.numeric(denominator.w)
mvp.w
## [,1]
## EEM -0.08352460
## EFA 0.01952054
## GLD 0.07442507
## IWM -0.15410052
## IYR -0.10224186
## QQQ -0.20279976
## SPY 1.03898880
## TLT 0.40973232
Q4
ret.sample.m <- ret.sample1 %>% to.monthly(indexAt = "lastof", OHLC = FALSE)
## Warning in to.period(x, "months", indexAt = indexAt, name = name, ...):
## missing values removed from data
avg.ret.m <- colMeans(coredata(ret.sample.m))
cov.d.m <- cov(coredata(ret.sample.m))
numerator.m <- solve(cov.d.m) %*% ones.vec
denominator.m <- t(ones.vec)%*%numerator.m
mvp.m <- numerator.m / as.numeric(denominator.m)
mvp.m
## [,1]
## EEM -0.13936904
## EFA -0.08733580
## GLD -0.01354903
## IWM -0.34925503
## IYR -0.04351678
## QQQ 0.04773161
## SPY 1.15577650
## TLT 0.42951757
Q5.
n <- length(tickers)
constraints = new.constraints(n, lb=-Inf, ub=+Inf)
constraints <- add.constraints(rep(1,n),1, type = "=", constraints)
ia.weekly <- create.historical.ia(ret.sample.w, 250)
weight.weekly <- min.risk.portfolio(ia.weekly,constraints)
## Loading required package: kernlab
##
## Attaching package: 'kernlab'
## The following object is masked _by_ '.GlobalEnv':
##
## cross
## The following object is masked from 'package:purrr':
##
## cross
## The following object is masked from 'package:ggplot2':
##
## alpha
ia.monthly <- create.historical.ia(ret.sample.m, 250)
weight.monthly <- min.risk.portfolio(ia.monthly,constraints)
Q6.
options(digits=5)
firm_data1 = read.csv("~/FamaFrench_mon_69_98_3stocks.csv")
head(firm_data1)
## date Mkt.RF SMB HML RF ge ibm mobil CRSP
## 1 196901 -1.20 -0.80 1.57 0.53 -1.1984 -5.9524 -1.4043 -0.6714
## 2 196902 -5.82 -3.90 0.93 0.46 -6.0377 -0.7004 -7.8431 -5.3641
## 3 196903 2.59 -0.28 -0.45 0.46 6.6474 7.0303 21.5130 3.0505
## 4 196904 1.52 -0.85 0.06 0.53 5.9621 4.4586 2.9961 2.0528
## 5 196905 0.02 -0.27 0.74 0.48 -3.5806 -2.5000 2.6667 0.5038
## 6 196906 -7.25 -5.31 -1.15 0.51 -3.8196 5.8777 -12.9870 -6.7388
attach(firm_data1)
library("Ecdat")
## Loading required package: Ecfun
##
## Attaching package: 'Ecfun'
## The following object is masked from 'package:base':
##
## sign
##
## Attaching package: 'Ecdat'
## The following object is masked from 'package:datasets':
##
## Orange
library("robust")
## Loading required package: fit.models
stocks5=cbind(N=Mkt.RF,SMB,HML)
fit = lm(cbind(Mkt.RF,SMB,HML)~RF)
options(digits=3)
pairs(cbind(Mkt.RF, SMB, HML))
cor(fit$residuals)
## Mkt.RF SMB HML
## Mkt.RF 1.000 0.323 -0.411
## SMB 0.323 1.000 -0.140
## HML -0.411 -0.140 1.000
covRob(fit$residuals,cor=F)
## Call:
## covRob(data = fit$residuals, corr = F)
##
## Robust Estimate of Covariance:
## Mkt.RF SMB HML
## Mkt.RF 16.02 3.48 -3.99
## SMB 3.48 6.08 -1.96
## HML -3.99 -1.96 5.86
##
## Robust Estimate of Location:
## Mkt.RF SMB HML
## 0.1124 0.0491 -0.2154
cor.test(fit$residuals[,1], fit$residuals[,2])
##
## Pearson's product-moment correlation
##
## data: fit$residuals[, 1] and fit$residuals[, 2]
## t = 6, df = 358, p-value = 3e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.227 0.413
## sample estimates:
## cor
## 0.323
cor.test(fit$residuals[,1], fit$residuals[,3])
##
## Pearson's product-moment correlation
##
## data: fit$residuals[, 1] and fit$residuals[, 3]
## t = -9, df = 358, p-value = 4e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.494 -0.322
## sample estimates:
## cor
## -0.411
cor.test(fit$residuals[,2], fit$residuals[,3])
##
## Pearson's product-moment correlation
##
## data: fit$residuals[, 2] and fit$residuals[, 3]
## t = -3, df = 358, p-value = 0.008
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2399 -0.0372
## sample estimates:
## cor
## -0.14
pairs(fit$residuals)
n=dim(firm_data1)[1]
sigF = as.matrix(var(cbind(Mkt.RF,SMB,HML)))
sigF
## Mkt.RF SMB HML
## Mkt.RF 21.15 4.23 -5.10
## SMB 4.23 8.18 -1.08
## HML -5.10 -1.08 7.18
bbeta = as.matrix(fit$coef)
bbeta = t( bbeta[-1,])
bbeta
## Mkt.RF SMB HML
## [1,] -2.66 -0.14 0.655
resig2 = apply((fit$resid)^2, 2, sum)/(n-3-1)
resig2 = diag(resig2)
cov_ff3 = sigF
cov_ff3
## Mkt.RF SMB HML
## Mkt.RF 21.15 4.23 -5.10
## SMB 4.23 8.18 -1.08
## HML -5.10 -1.08 7.18
cov2cor(cov_ff3)
## Mkt.RF SMB HML
## Mkt.RF 1.000 0.322 -0.414
## SMB 0.322 1.000 -0.140
## HML -0.414 -0.140 1.000
cov_hist = cov(stocks5)
cov2cor(cov_hist)
## N SMB HML
## N 1.000 0.322 -0.414
## SMB 0.322 1.000 -0.140
## HML -0.414 -0.140 1.000
one.vec=rep(1,3)
a = solve(cov_ff3)%*%one.vec
b = t(one.vec)%*%a
mvp.w =a / as.numeric(b)
mvp.w
## [,1]
## Mkt.RF 0.191
## SMB 0.279
## HML 0.529
a1 = solve(cov_hist)%*%one.vec
b1 = t(one.vec)%*%a1
mvp.w1 =a1 / as.numeric(b1)
mvp.w1
## [,1]
## N 0.191
## SMB 0.279
## HML 0.529