library(tidyverse)
## -- Attaching packages -------------------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1 v purrr 0.3.3
## v tibble 2.1.3 v dplyr 0.8.3
## v tidyr 1.0.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.4.0
## -- Conflicts ----------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(xts)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo
##
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
##
## first, last
library(PerformanceAnalytics)
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
## The following object is masked from 'package:purrr':
##
## compact
library(quantmod)
## 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(fPortfolio)
## Warning: package 'fPortfolio' was built under R version 3.6.3
## Loading required package: timeDate
##
## Attaching package: 'timeDate'
## The following objects are masked from 'package:PerformanceAnalytics':
##
## kurtosis, skewness
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 3.6.3
##
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
##
## time<-
## Loading required package: fBasics
## Warning: package 'fBasics' was built under R version 3.6.3
##
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
##
## volatility
## Loading required package: fAssets
## Warning: package 'fAssets' was built under R version 3.6.3
library(readr)
library(fit.models)
## Warning: package 'fit.models' was built under R version 3.6.3
library(robust)
## Warning: package 'robust' was built under R version 3.6.3
library(Ecdat)
## Warning: package 'Ecdat' was built under R version 3.6.3
## Loading required package: Ecfun
## Warning: package 'Ecfun' was built under R version 3.6.3
##
## Attaching package: 'Ecfun'
## The following object is masked from 'package:base':
##
## sign
##
## Attaching package: 'Ecdat'
## The following object is masked from 'package:datasets':
##
## Orange
tickers <- c("SPY", "QQQ", "EEM","IWM","EFA","TLT","IYR","GLD")
data <- new.env()
getSymbols(tickers, src = 'yahoo', from = "2010-01-01", to = "2018-12-31", 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"
#Question 1
portfolioPrices <- NULL
for(ticker in tickers) {
portfolioPrices <- cbind(portfolioPrices,
getSymbols(ticker, src = "yahoo", from='2010-01-01',to = "2018-12-31", periodicity = 'daily', auto.assign=FALSE)[,4])
}
portfolioPrices <- portfolioPrices[apply(portfolioPrices,1,function(x) all(!is.na(x))),]
colnames(portfolioPrices) <- tickers
portfolioReturns <- na.omit(ROC(portfolioPrices, type="discrete"))
portfolioReturns <- as.timeSeries(portfolioReturns)
colnames(portfolioReturns) <- tickers
portfolioReturns <- as.timeSeries(portfolioReturns)
head(portfolioReturns)
## GMT
## SPY QQQ EEM IWM
## 2010-01-05 0.002647093 0.0000000000 0.007258277 -0.0034386058
## 2010-01-06 0.000704057 -0.0060318615 0.002092073 -0.0009409818
## 2010-01-07 0.004221291 0.0006501734 -0.005799118 0.0073782887
## 2010-01-08 0.003327769 0.0082304747 0.007932804 0.0054542467
## 2010-01-11 0.001396552 -0.0040815898 -0.002083333 -0.0040296809
## 2010-01-12 -0.009326235 -0.0125108280 -0.016005636 -0.0108932614
## EFA TLT IYR GLD
## 2010-01-05 0.0008813503 0.0064580894 0.002401157 -0.0009108014
## 2010-01-06 0.0042268581 -0.0133864256 -0.000435453 0.0164995902
## 2010-01-07 -0.0038583129 0.0016820139 0.008932440 -0.0061878037
## 2010-01-08 0.0079225530 -0.0004477891 -0.006694040 0.0049630301
## 2010-01-11 0.0082096245 -0.0054877141 0.004782630 0.0132889913
## 2010-01-12 -0.0117810116 0.0171170828 -0.016875854 -0.0209127164
#Question 2
or_matrix <- cor(portfolioReturns)
cov_matrix <- cov(portfolioReturns)
write.csv(cov_matrix, "covmatrix.csv")
mvp <- minvariancePortfolio(portfolioReturns, spec = portfolioSpec(), constraints = "Longonly")
mvpweights <- getWeights(mvp)
mvpweights
## SPY QQQ EEM IWM EFA TLT IYR
## 0.4259534 0.0000000 0.0000000 0.0000000 0.0000000 0.4495901 0.0000000
## GLD
## 0.1244564
#Question 3
portfolioPrices <- NULL
for(ticker in tickers) {
portfolioPrices <- cbind(portfolioPrices,
getSymbols.yahoo(ticker, from='2010-01-01',to = "2013-12-31", periodicity = 'weekly', auto.assign=FALSE)[,4])
}
portfolioPrices <- portfolioPrices[apply(portfolioPrices,1,function(x) all(!is.na(x))),]
colnames(portfolioPrices) <- tickers
portfolioReturns <- na.omit(ROC(portfolioPrices, type="discrete"))
portfolioReturns <- as.timeSeries(portfolioReturns)
colnames(portfolioReturns) <- tickers
portfolioReturns <- as.timeSeries(portfolioReturns)
head(portfolioReturns)
## GMT
## SPY QQQ EEM IWM EFA
## 2010-01-08 0.00648041 0.004765021 -0.009799393 0.007012701 0.02042254
## 2010-01-15 -0.02810409 -0.019400669 -0.045475991 -0.027236179 -0.05210490
## 2010-01-22 -0.02802146 -0.042646800 -0.044680254 -0.033089436 -0.03076081
## 2010-01-29 -0.01961866 -0.021354765 -0.027907028 -0.028627806 -0.03154930
## 2010-02-05 0.01587744 0.024636298 0.034821931 0.026084027 0.01221642
## 2010-02-12 0.02570986 0.027020839 0.020292859 0.039617001 0.02049806
## TLT IYR GLD
## 2010-01-08 0.013321336 -0.003886849 0.010918598
## 2010-01-15 0.018117621 -0.028181205 -0.041595966
## 2010-01-22 -0.007052973 -0.024760226 -0.008289094
## 2010-01-29 0.002950464 -0.021729208 -0.019815927
## 2010-02-05 -0.021355404 0.001870423 0.026444322
## 2010-02-12 -0.009463360 0.059976666 0.026603249
or_matrix <- cor(portfolioReturns)
cov_matrix <- cov(portfolioReturns)
write.csv(cov_matrix, "covmatrix.csv")
mvp <- minvariancePortfolio(portfolioReturns, spec = portfolioSpec(), constraints = "Longonly")
mvpweights <- getWeights(mvp)
mvpweights
## SPY QQQ EEM IWM EFA TLT
## 0.49303564 0.00000000 0.00000000 0.00000000 0.00000000 0.48111726
## IYR GLD
## 0.00000000 0.02584709
#Question 4
portfolioPrices <- NULL
for(ticker in tickers) {
portfolioPrices <- cbind(portfolioPrices,
getSymbols.yahoo(ticker, from='2010-01-01',to = "2013-12-31", periodicity = 'monthly', auto.assign=FALSE)[,4])
}
portfolioPrices <- portfolioPrices[apply(portfolioPrices,1,function(x) all(!is.na(x))),]
colnames(portfolioPrices) <- tickers
portfolioReturns <- na.omit(ROC(portfolioPrices, type="discrete"))
portfolioReturns <- as.timeSeries(portfolioReturns)
colnames(portfolioReturns) <- tickers
portfolioReturns <- as.timeSeries(portfolioReturns)
head(portfolioReturns)
## GMT
## SPY QQQ EEM IWM EFA
## 2010-02-01 0.03119470 0.04603872 0.017763846 0.04475126 0.002667664
## 2010-03-01 0.05652883 0.07596073 0.081108832 0.07961790 0.063854068
## 2010-04-01 0.01547007 0.02242529 -0.001661918 0.05678464 -0.028045731
## 2010-05-01 -0.07945455 -0.07392372 -0.093935817 -0.07536639 -0.111927954
## 2010-06-01 -0.05623116 -0.06337717 -0.020472390 -0.07743398 -0.037458651
## 2010-07-01 0.06830068 0.07258258 0.109324812 0.06380887 0.116104112
## TLT IYR GLD
## 2010-02-01 -0.00693316 0.05457055 0.032748217
## 2010-03-01 -0.02367185 0.08689957 -0.004386393
## 2010-04-01 0.02938544 0.06388108 0.058834366
## 2010-05-01 0.04743301 -0.05683531 0.030513141
## 2010-06-01 0.05440415 -0.05485489 0.023553189
## 2010-07-01 -0.01248154 0.09404794 -0.050871154
or_matrix <- cor(portfolioReturns)
cov_matrix <- cov(portfolioReturns)
write.csv(cov_matrix, "covmatrix.csv")
mvp <- minvariancePortfolio(portfolioReturns, spec = portfolioSpec(), constraints = "Longonly")
mvpweights <- getWeights(mvp)
mvpweights
## SPY QQQ EEM IWM EFA TLT
## 0.44235680 0.00000000 0.00000000 0.04949215 0.00000000 0.48955621
## IYR GLD
## 0.00000000 0.01859485
retdata = read.csv("FFResearchDataFactors.CSV")
attach(retdata)
stocks5=cbind(N=Mkt.RF,SMB,HML)
fit = lm(cbind(Mkt.RF,SMB,HML)~RF)
options(digits=3)
head(SMB)
## [1] -2.3 -1.4 -1.32 0.04 -0.2 -0.04
## 758 Levels: -0.01 -0.03 -0.04 -0.05 -0.07 -0.08 -0.09 -0.1 -0.11 ... SMB
pairs(cbind(Mkt.RF,SMB,HML))
cor(fit$residuals)
## Mkt.RF SMB HML
## Mkt.RF 1.00000 0.1982 0.00895
## SMB 0.19824 1.0000 0.08089
## HML 0.00895 0.0809 1.00000
covRob(fit$residuals,cor=F)
## Call:
## covRob(data = fit$residuals, corr = F)
##
## Robust Estimate of Covariance:
## Mkt.RF SMB HML
## Mkt.RF 66724 15796 -1920
## SMB 15796 48106 2864
## HML -1920 2864 44271
##
## Robust Estimate of Location:
## Mkt.RF SMB HML
## -0.0452 -3.4413 -2.6780
cor.test(fit$residuals[,1], fit$residuals[,2])
##
## Pearson's product-moment correlation
##
## data: fit$residuals[, 1] and fit$residuals[, 2]
## t = 7, df = 1224, p-value = 2e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.144 0.251
## sample estimates:
## cor
## 0.198
cor.test(fit$residuals[,1], fit$residuals[,3])
##
## Pearson's product-moment correlation
##
## data: fit$residuals[, 1] and fit$residuals[, 3]
## t = 0.3, df = 1224, p-value = 0.8
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0471 0.0649
## sample estimates:
## cor
## 0.00895
cor.test(fit$residuals[,2], fit$residuals[,3])
##
## Pearson's product-moment correlation
##
## data: fit$residuals[, 2] and fit$residuals[, 3]
## t = 3, df = 1224, p-value = 0.005
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.025 0.136
## sample estimates:
## cor
## 0.0809
pairs(fit$residuals)
n=dim(retdata)[1]
sigF = as.matrix(var(cbind(Mkt.RF,SMB,HML)))
sigF
## Mkt.RF SMB HML
## Mkt.RF 65970 12827 1158
## SMB 12827 47529 4790
## HML 1158 4790 46390
bbeta = as.matrix(fit$coef)
bbeta = t( bbeta[-1,])
bbeta
## RF-0.01 RF-0.02 RF-0.03 RF-0.06 RF0 RF0.01 RF0.02 RF0.03 RF0.04
## Mkt.RF 333 417 135 95 435 505 448 496 507
## SMB 276 545 218 207 339 414 422 351 320
## HML 407 212 214 105 229 374 376 376 338
## RF0.05 RF0.06 RF0.07 RF0.08 RF0.09 RF0.1 RF0.11 RF0.12 RF0.13
## Mkt.RF 513 365 484 513 498 538 500 321 489
## SMB 475 280 424 237 339 320 331 346 321
## HML 401 335 363 347 262 341 311 209 344
## RF0.14 RF0.15 RF0.16 RF0.17 RF0.18 RF0.19 RF0.2 RF0.21 RF0.22
## Mkt.RF 403 489 367 444 489 422 478 429 472
## SMB 386 356 270 365 346 255 265 299 345
## HML 271 355 298 289 235 347 278 268 304
## RF0.23 RF0.24 RF0.25 RF0.26 RF0.27 RF0.28 RF0.29 RF0.3 RF0.31
## Mkt.RF 440 361 433 366 336 488 468 509 428
## SMB 401 291 279 272 310 336 249 314 306
## HML 354 324 341 260 361 387 268 386 316
## RF0.32 RF0.33 RF0.34 RF0.35 RF0.36 RF0.37 RF0.38 RF0.39 RF0.4
## Mkt.RF 489 491 278 333 340 363 362 488 395
## SMB 330 502 356 447 397 353 304 326 342
## HML 314 339 336 398 398 271 332 280 358
## RF0.41 RF0.42 RF0.43 RF0.44 RF0.45 RF0.46 RF0.47 RF0.48 RF0.49
## Mkt.RF 496 433 541 505 510 452 501 332 523
## SMB 402 331 343 334 428 337 311 253 413
## HML 265 314 313 259 417 353 397 443 266
## RF0.5 RF0.51 RF0.52 RF0.53 RF0.54 RF0.55 RF0.56 RF0.57 RF0.58
## Mkt.RF 540 260 482 452 537 445 455 368 329
## SMB 375 299 286 297 442 202 432 254 372
## HML 493 402 368 403 292 440 321 431 433
## RF0.59 RF0.6 RF0.61 RF0.62 RF0.63 RF0.64 RF0.65 RF0.66 RF0.67
## Mkt.RF 369 423 508 341 381 478 475 446 569
## SMB 462 283 378 302 590 422 240 223 233
## HML 397 388 314 376 338 273 289 274 278
## RF0.68 RF0.69 RF0.7 RF0.71 RF0.72 RF0.73 RF0.74 RF0.75 RF0.76
## Mkt.RF 366 438 471 396 49 209 361 404 261
## SMB 331 365 393 209 372 211 444 191 259
## HML 285 123 322 520 635 489 551 278 597
## RF0.77 RF0.78 RF0.79 RF0.8 RF0.81 RF0.82 RF0.83 RF0.86 RF0.87
## Mkt.RF 654 370 669 410 481 169 278 408 287
## SMB 554 414 3 367 412 277 22 277 356
## HML 389 276 73 425 445 418 113 644 279
## RF0.89 RF0.92 RF0.95 RF0.96 RF0.98 RF0.99 RF1 RF1.02 RF1.04 RF1.05
## Mkt.RF 63 302 470 487 101 789 44 708 280 211
## SMB 144 404 621 329 16 591 98 614 615 433
## HML 397 710 186 436 636 242 387 692 721 356
## RF1.06 RF1.07 RF1.08 RF1.1 RF1.13 RF1.15 RF1.2 RF1.21 RF1.24 RF1.26
## Mkt.RF 240 453 150 550 664 379 577 526 208 706
## SMB 402 250 684 661 489 536 569 509 196 446
## HML 493 366 552 291 218 37 667 257 372 434
## RF1.28 RF1.31 RF1.35 RF1.49 RF1.54 RF1.57 RF1.6 RF1.65 RF1.66
## Mkt.RF 332 264 162 551 774 627 244 199 529
## SMB 154 25 72 312 529 338 626 691 337
## HML 675 576 684 312 517 703 412 729 617
## RF1.81 RF1.82 RF10.38 RF10.54 RF11.24 RF14.71 RF2.13 RF2.14 RF2.41
## Mkt.RF 328 52 536 517 626 143 630 635 241
## SMB 261 95 607 751 714 739 487 330 318
## HML 336 325 184 518 230 597 686 165 167
## RF2.46 RF2.66 RF2.73 RF2.9 RF2.95 RF2.98 RF3.12 RF3.14 RF3.51
## Mkt.RF 815 78 126 879 894 649 594 128 831
## SMB 96 219 349 718 712 152 263 216 743
## HML 133 297 742 529 497 738 392 318 595
## RF3.54 RF3.56 RF3.83 RF3.84 RF3.9 RF3.93 RF4.21 RF4.39 RF4.66
## Mkt.RF 533 714 136 535 250 515 629 528 444
## SMB 82 683 535 159 88 609 724 715 344
## HML 750 310 531 487 61 726 332 163 173
## RF4.68 RF4.75 RF4.76 RF4.8 RF4.86 RF5.08 RF5.12 RF5.21 RF5.26
## Mkt.RF 621 144 132 516 552 625 353 713 632
## SMB 530 281 580 380 234 528 610 439 339
## HML 270 515 60 520 335 596 730 619 526
## RF5.47 RF5.6 RF5.8 RF5.89 RF6.16 RF6.35 RF6.52 RF6.58 RF6.93 RF7.18
## Mkt.RF 234 673 711 141 510 525 321 140 202 443
## SMB 157 442 531 122 356 721 158 161 233 527
## HML 136 430 744 648 745 521 593 340 527 375
## RF7.72 RF7.81 RF8 RF8.37 RF8.8 RF9.85 RFRF
## Mkt.RF 631 134 243 620 540 311 895
## SMB 370 163 52 160 523 351 757
## HML 455 339 746 285 591 532 751
resig2 = apply((fit$resid)^2, 2, sum)/(n-3-1)
resig2 = diag(resig2)
cov_ff3 = sigF
cov_ff3
## Mkt.RF SMB HML
## Mkt.RF 65970 12827 1158
## SMB 12827 47529 4790
## HML 1158 4790 46390
cov2cor(cov_ff3)
## Mkt.RF SMB HML
## Mkt.RF 1.0000 0.229 0.0209
## SMB 0.2291 1.000 0.1020
## HML 0.0209 0.102 1.0000
cov_hist = cov(stocks5)
cov2cor(cov_hist)
## N SMB HML
## N 1.0000 0.229 0.0209
## SMB 0.2291 1.000 0.1020
## HML 0.0209 0.102 1.0000
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.248
## SMB 0.337
## HML 0.415
a1 = solve(cov_hist)%*%one.vec
b1 = t(one.vec)%*%a1
mvp.w1 =a1 / as.numeric(b1)
mvp.w1
## [,1]
## N 0.248
## SMB 0.337
## HML 0.415