rm(list = ls())
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(timetk)
library(purrr)
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
library(lpSolve)
library(SIT)
## Loading required package: SIT.date
##
## Attaching package: 'SIT'
## The following object is masked from 'package:TTR':
##
## DVI
## The following object is masked from 'package:purrr':
##
## cross
## The following object is masked from 'package:base':
##
## close
library(PerformanceAnalytics)
##
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
##
## legend
library(xts)
library(tidyr)
Selected stocks
- HollyFrontier Corporation (HFC)
- Kansas City Southern (KSU)
- Jack Henry & Associates, Inc. (JKHY)
- Cerner Corporation (CERN)
- NetEase, Inc. (NTES)
- NVIDIA Corporation (NVDA)
- SBA Communications Corporation (SBAC)
- Equinix, Inc. (EQIX)
Import data from Yahoo Finance
tickers = spl('HFC,KSU,JKHY,CERN,NTES,NVDA,SBAC,EQIX')
datas <- new.env()
getSymbols(tickers, src = 'yahoo', from = '1980-01-01', env = datas, 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] "HFC" "KSU" "JKHY" "CERN" "NTES" "NVDA" "SBAC" "EQIX"
for(i in ls(datas)) datas[[i]] = adjustOHLC(datas[[i]], use.Adjusted=T)
names(datas)
## [1] "EQIX" "CERN" "NTES" "HFC" ".getSymbols"
## [6] "JKHY" "NVDA" "KSU" "SBAC"
head(datas$EQIX)
## EQIX.Open EQIX.High EQIX.Low EQIX.Close EQIX.Volume EQIX.Adjusted
## 2000-08-11 331.0727 363.5308 311.5979 340.8102 792528 340.8102
## 2000-08-14 339.9987 340.8102 311.5979 311.5979 80228 311.5979
## 2000-08-15 315.6551 321.3353 311.5979 314.8437 48422 314.8437
## 2000-08-16 313.2208 339.1873 313.2208 321.3353 42166 321.3353
## 2000-08-17 318.0895 324.5811 309.9750 311.5979 47247 311.5979
## 2000-08-18 325.3925 373.2682 321.3353 371.6453 60400 371.6453
tail(datas$EQIX)
## EQIX.Open EQIX.High EQIX.Low EQIX.Close EQIX.Volume EQIX.Adjusted
## 2021-06-15 819.76 826.280 807.03 809.00 526300 809.00
## 2021-06-16 813.00 819.680 806.75 809.97 636400 809.97
## 2021-06-17 809.97 816.820 805.90 815.46 661400 815.46
## 2021-06-18 815.50 829.610 799.84 822.77 3751600 822.77
## 2021-06-21 818.67 826.810 802.65 822.50 1011300 822.50
## 2021-06-22 823.65 824.985 811.57 815.72 207614 815.72
bt.prep(datas, align='remove.na')
names(datas)
## [1] "EQIX" "prices" "CERN" "NTES"
## [5] "dates" "weight" "HFC" ".getSymbols"
## [9] "symbolnames" "execution.price" "JKHY" "NVDA"
## [13] "KSU" "SBAC"
head(datas$NTES)
## NTES.Open NTES.High NTES.Low NTES.Close NTES.Volume
## 2000-08-11 08:00:00 0.2552053 0.2552053 0.2229689 0.231028 984000
## 2000-08-14 08:00:00 0.2364010 0.2390874 0.2095372 0.214910 12756000
## 2000-08-15 08:00:00 0.2135672 0.2189400 0.2122236 0.214910 1364000
## 2000-08-16 08:00:00 0.2149100 0.2256555 0.2122236 0.214910 2280000
## 2000-08-17 08:00:00 0.2256557 0.2498331 0.2122238 0.241774 1838000
## 2000-08-18 08:00:00 0.2310285 0.2632650 0.2256557 0.263265 2378000
## NTES.Adjusted
## 2000-08-11 08:00:00 0.231028
## 2000-08-14 08:00:00 0.214910
## 2000-08-15 08:00:00 0.214910
## 2000-08-16 08:00:00 0.214910
## 2000-08-17 08:00:00 0.241774
## 2000-08-18 08:00:00 0.263265
head(datas$prices)
## CERN EQIX HFC JKHY KSU NTES
## 2000-08-11 08:00:00 4.534172 340.8102 0.380951 16.14683 5.983018 0.231028
## 2000-08-14 08:00:00 4.602987 311.5979 0.382915 16.41595 5.983018 0.214910
## 2000-08-15 08:00:00 4.488295 314.8437 0.388806 16.68506 5.983018 0.214910
## 2000-08-16 08:00:00 4.610633 321.3353 0.392733 16.68506 6.038935 0.214910
## 2000-08-17 08:00:00 4.549465 311.5979 0.396660 16.85631 7.157254 0.241774
## 2000-08-18 08:00:00 4.450064 371.6453 0.412370 16.63613 7.548665 0.263265
## NVDA SBAC
## 2000-08-11 08:00:00 9.868440 44.39121
## 2000-08-14 08:00:00 9.964158 45.06941
## 2000-08-15 08:00:00 10.904580 44.26790
## 2000-08-16 08:00:00 11.179769 44.14460
## 2000-08-17 08:00:00 11.390343 44.82280
## 2000-08-18 08:00:00 11.792356 44.51452
Convert prices into monthly data
prices_monthly <- to.monthly(datas$prices, indexAt = "last", OHLC = FALSE) # indexAt = 'endof'
head(prices_monthly)
## CERN EQIX HFC JKHY KSU NTES NVDA
## 2000-08-31 4.656511 420.3326 0.404515 17.37007 8.275572 0.327738 12.156082
## 2000-09-29 5.681096 230.4527 0.394430 16.99795 7.772331 0.233715 12.538951
## 2000-10-31 7.577343 253.1732 0.398413 21.55360 7.772331 0.150437 9.516680
## 2000-11-30 5.979296 113.6034 0.426302 21.03472 7.548665 0.150437 6.202473
## 2000-12-29 5.658157 113.6034 0.603876 24.36888 9.058400 0.131632 5.017973
## 2001-01-31 5.910480 154.1760 0.579721 17.89667 11.299516 0.096709 7.906240
## SBAC
## 2000-08-31 44.02129
## 2000-09-29 41.37014
## 2000-10-31 49.44688
## 2000-11-30 36.74606
## 2000-12-29 40.50699
## 2001-01-31 44.69948
tail(prices_monthly)
## CERN EQIX HFC JKHY KSU NTES NVDA
## 2021-01-29 79.86831 733.8204 28.20067 143.9136 201.7551 114.6340 519.2917
## 2021-02-26 68.93140 645.7377 37.88000 147.5415 211.3815 109.5099 548.2651
## 2021-03-31 71.66313 676.8622 35.78000 151.2740 263.4353 102.9988 533.8077
## 2021-04-30 75.05000 717.8670 35.00000 162.3514 291.6733 111.7765 600.2424
## 2021-05-28 78.25000 736.7200 32.47000 154.1500 297.1333 117.6316 649.6312
## 2021-06-22 78.95000 815.7200 33.84000 166.5800 284.0100 107.2000 753.6800
## SBAC
## 2021-01-29 267.4947
## 2021-02-26 254.0139
## 2021-03-31 276.9974
## 2021-04-30 299.1233
## 2021-05-28 298.1200
## 2021-06-22 320.5800
Calculate monthly returns
monthly.ret <- na.omit(Return.calculate(prices_monthly, method = "discrete"))
head(monthly.ret)
## CERN EQIX HFC JKHY KSU
## 2000-09-29 0.22003277 -0.45173732 -0.02493109 -0.021422832 -0.06081042
## 2000-10-31 0.33378190 0.09859117 0.01009812 0.268011437 0.00000000
## 2000-11-30 -0.21089807 -0.55128194 0.07000023 -0.024073896 -0.02877721
## 2000-12-29 -0.05370850 0.00000000 0.41654508 0.158507916 0.20000026
## 2001-01-31 0.04459456 0.35714246 -0.03999993 -0.265593563 0.24740749
## 2001-02-28 0.05950853 -0.36842076 0.15833306 0.002605513 0.18606487
## NTES NVDA SBAC
## 2000-09-29 -0.2868846 0.03149609 -0.06022409
## 2000-10-31 -0.3563229 -0.24103061 0.19523103
## 2000-11-30 0.0000000 -0.34825244 -0.25685788
## 2000-12-29 -0.1250025 -0.19097221 0.10234910
## 2001-01-31 -0.2653078 0.57558440 0.10350059
## 2001-02-28 -0.1388806 -0.13438259 -0.27172404
tail(monthly.ret)
## CERN EQIX HFC JKHY KSU
## 2021-01-29 0.020769645 0.03609735 0.10096708 -0.10617956 -0.007152258
## 2021-02-26 -0.136936729 -0.12003306 0.34323077 0.02520909 0.047713000
## 2021-03-31 0.039629658 0.04820003 -0.05543828 0.02529826 0.246255333
## 2021-04-30 0.047260996 0.06058066 -0.02179986 0.07322702 0.107191544
## 2021-05-28 0.042638199 0.02626248 -0.07228569 -0.05051636 0.018719438
## 2021-06-22 0.008945649 0.10723206 0.04219276 0.08063580 -0.044166342
## NTES NVDA SBAC
## 2021-01-29 0.20068914 -0.004998005 -0.047708523
## 2021-02-26 -0.04469953 0.055793931 -0.050396423
## 2021-03-31 -0.05945701 -0.026369358 0.090481479
## 2021-04-30 0.08522173 0.124454474 0.079877468
## 2021-05-28 0.05238264 0.082281309 -0.003354122
## 2021-06-22 -0.08868063 0.160166004 0.075338764
download Fama-French factors from website
url <- "http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/ftp/F-F_Research_Data_Factors_daily_CSV.zip"
temp <- tempfile()
download.file(url, temp, method = "libcurl", mode = "wb")
unzip(temp, "F-F_Research_Data_Factors_daily.CSV")
unlink(temp)
mydata <- read.csv("F-F_Research_Data_Factors_daily.CSV", skip = 4)
mydata <- mydata[-nrow(mydata), ] # remove last row
fama_lib <- xts(x = mydata[, c(2,3,4)], order.by = as.Date(paste(mydata[, 1]), "%Y%m%d"))
str(fama_lib)
## An 'xts' object on 1926-07-01/2021-04-30 containing:
## Data: num [1:24978, 1:3] 0.1 0.45 0.17 0.09 0.21 -0.71 0.62 0.04 0.48 0.04 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:3] "Mkt.RF" "SMB" "HML"
## Indexed by objects of class: [Date] TZ: UTC
## xts Attributes:
## NULL
head(fama_lib)
## Mkt.RF SMB HML
## 1926-07-01 0.10 -0.24 -0.28
## 1926-07-02 0.45 -0.32 -0.08
## 1926-07-06 0.17 0.27 -0.35
## 1926-07-07 0.09 -0.59 0.03
## 1926-07-08 0.21 -0.36 0.15
## 1926-07-09 -0.71 0.44 0.56
dates <- '2005-01::2021-03'
X <- monthly.ret[dates]
head(X)
## CERN EQIX HFC JKHY KSU
## 2005-01-31 -0.064321874 -0.01871805 0.085755326 0.04419880 -0.015228971
## 2005-02-28 0.046834207 0.04411059 0.245870052 -0.04304979 0.123711670
## 2005-03-31 0.008256439 -0.03311266 -0.009247788 -0.09370291 -0.018348544
## 2005-04-29 0.105503646 -0.17383080 -0.080493556 -0.04446934 -0.017653324
## 2005-05-31 0.125753635 0.08919395 0.115553438 0.03170599 0.056553937
## 2005-06-30 0.040092075 0.13753285 0.223531525 0.03504789 0.009504934
## NTES NVDA SBAC
## 2005-01-31 -0.20351472 -0.02716471 -0.07758616
## 2005-02-28 0.00711766 0.26483441 0.01635514
## 2005-03-31 0.13568941 -0.18040737 0.05057478
## 2005-04-29 0.02447587 -0.07702011 -0.07221048
## 2005-05-31 0.04555602 0.23575029 0.32547203
## 2005-06-30 0.10592582 -0.01402225 0.20106786
tail(X)
## CERN EQIX HFC JKHY KSU
## 2020-10-30 -0.03043305 -0.03800675 -0.06088281 -0.088197332 -0.025935920
## 2020-11-30 0.06777003 -0.04237073 0.28327186 0.085059133 0.056943362
## 2020-12-31 0.05161503 0.02348847 0.10517315 0.009723934 0.098891057
## 2021-01-29 0.02076964 0.03609735 0.10096708 -0.106179561 -0.007152258
## 2021-02-26 -0.13693673 -0.12003306 0.34323077 0.025209088 0.047713000
## 2021-03-31 0.03962966 0.04820003 -0.05543828 0.025298265 0.246255333
## NTES NVDA SBAC
## 2020-10-30 -0.04557146 -0.073648422 -0.088263034
## 2020-11-30 0.04124882 0.069211734 -0.009444331
## 2020-12-31 0.06210449 -0.025567475 -0.017584742
## 2021-01-29 0.20068914 -0.004998005 -0.047708523
## 2021-02-26 -0.04469953 0.055793931 -0.050396423
## 2021-03-31 -0.05945701 -0.026369358 0.090481479
dim(X)
## [1] 195 8
fama_lib_month <- to.monthly(fama_lib, indexAt = "last", OHLC = FALSE) # indexAt = 'endof'
head(fama_lib_month)
## Mkt.RF SMB HML
## 1926-07-31 0.46 -0.12 -0.17
## 1926-08-31 0.58 -0.52 -0.04
## 1926-09-30 0.30 -0.22 0.06
## 1926-10-30 -0.25 0.20 0.25
## 1926-11-30 0.26 0.03 -0.50
## 1926-12-31 0.33 0.30 0.67
f <- fama_lib_month[dates]/100
head(f)
## Mkt.RF SMB HML
## 2005-01-31 0.0097 0.0076 0.0014
## 2005-02-28 -0.0061 0.0027 0.0007
## 2005-03-31 0.0000 -0.0006 0.0065
## 2005-04-29 0.0104 -0.0040 -0.0003
## 2005-05-31 -0.0049 0.0046 0.0033
## 2005-06-30 -0.0063 0.0021 0.0025
dim(f)
## [1] 195 3
======================================================================
Use Professor Palomar’s package to compute estimated coefficients
Based on CAPM model, compute covariance matrix for the 8-asset portfolio by using past 60-month returns from 2010/01 - 2014/12.
======================================================================
insample_range <- '2010-01::2014-12'
one_factor_model <- factorModel(X[insample_range], type = "M", econ_fact = f$Mkt.RF[insample_range])
names(one_factor_model)
## [1] "alpha" "beta" "factors" "residual"
cbind(alpha = one_factor_model$alpha, beta = one_factor_model$beta)
## alpha Mkt.RF
## CERN 0.021036176 -0.06045842
## EQIX 0.016959077 1.33127873
## HFC 0.028343752 0.51556901
## JKHY 0.018729683 0.11091914
## KSU 0.025434263 0.59261756
## NTES 0.021397145 -0.16540147
## NVDA 0.009871609 2.37834017
## SBAC 0.021622073 1.03437401
======================================================================
we can do the fitting using a more compact matrix notation
F_ <- cbind(ones = 1, f$Mkt.RF[insample_range])
Gamma <- t(X[insample_range]) %*% F_ %*% solve(t(F_) %*% F_) # better: Gamma <- t(solve(t(F_) %*% F_, t(F_) %*% X))
colnames(Gamma) <- c("alpha", "beta")
alpha <- Gamma[, 1] # or alpha <- Gamma[, "alpha"]
beta <- Gamma[, 2] # or beta <- Gamma[, "beta"]
print(Gamma)
## alpha beta
## CERN 0.021036176 -0.06045842
## EQIX 0.016959077 1.33127873
## HFC 0.028343752 0.51556901
## JKHY 0.018729683 0.11091914
## KSU 0.025434263 0.59261756
## NTES 0.021397145 -0.16540147
## NVDA 0.009871609 2.37834017
## SBAC 0.021622073 1.03437401
compute Sigma
X <- X[insample_range]
f <- f$Mkt.RF[insample_range]
T <- dim(X)[1]
E <- xts(t(t(X) - Gamma %*% t(F_)), index(X)) # residuals
Psi <- (1/(T-2)) * t(E) %*% E
Sigma <- as.numeric(var(f)) * beta %o% beta + diag(diag(Psi))
Sigma # This is covariance matrix computed from single factor model
## CERN EQIX HFC JKHY KSU
## CERN 3.897737e-03 -8.791295e-06 -3.404636e-06 -7.324709e-07 -3.913437e-06
## EQIX -8.791295e-06 5.703519e-03 7.496920e-05 1.612882e-05 8.617288e-05
## HFC -3.404636e-06 7.496920e-05 1.033759e-02 6.246265e-06 3.337247e-05
## JKHY -7.324709e-07 1.612882e-05 6.246265e-06 1.772552e-03 7.179730e-06
## KSU -3.913437e-06 8.617288e-05 3.337247e-05 7.179730e-06 5.689842e-03
## NTES 1.092253e-06 -2.405113e-05 -9.314366e-06 -2.003886e-06 -1.070634e-05
## NVDA -1.570572e-05 3.458359e-04 1.339331e-04 2.881427e-05 1.539485e-04
## SBAC -6.830641e-06 1.504090e-04 5.824941e-05 1.253173e-05 6.695441e-05
## NTES NVDA SBAC
## CERN 1.092253e-06 -1.570572e-05 -6.830641e-06
## EQIX -2.405113e-05 3.458359e-04 1.504090e-04
## HFC -9.314366e-06 1.339331e-04 5.824941e-05
## JKHY -2.003886e-06 2.881427e-05 1.253173e-05
## KSU -1.070634e-05 1.539485e-04 6.695441e-05
## NTES 8.817481e-03 -4.296754e-05 -1.868719e-05
## NVDA -4.296754e-05 1.441872e-02 2.687068e-04
## SBAC -1.868719e-05 2.687068e-04 2.493187e-03
=======================================================================
We can also use lm() to compute estimated coefficients
fit = lm(formula = X~f)
sigF = as.numeric(var(f))
beta_ = as.matrix(fit$coefficients)
beta_ = as.matrix(beta_[-1,])
beta_
## [,1]
## CERN -0.06045842
## EQIX 1.33127873
## HFC 0.51556901
## JKHY 0.11091914
## KSU 0.59261756
## NTES -0.16540147
## NVDA 2.37834017
## SBAC 1.03437401
sigeps = crossprod(fit$residuals)/(T-2)
# sigeps = as.matrix(var(fit$residuals)) # you can use this way too
sigeps = diag(diag(sigeps))
sigeps
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.003897338 0.000000000 0.00000000 0.000000000 0.000000000 0.000000000
## [2,] 0.000000000 0.005509937 0.00000000 0.000000000 0.000000000 0.000000000
## [3,] 0.000000000 0.000000000 0.01030856 0.000000000 0.000000000 0.000000000
## [4,] 0.000000000 0.000000000 0.00000000 0.001771208 0.000000000 0.000000000
## [5,] 0.000000000 0.000000000 0.00000000 0.000000000 0.005651482 0.000000000
## [6,] 0.000000000 0.000000000 0.00000000 0.000000000 0.000000000 0.008814493
## [7,] 0.000000000 0.000000000 0.00000000 0.000000000 0.000000000 0.000000000
## [8,] 0.000000000 0.000000000 0.00000000 0.000000000 0.000000000 0.000000000
## [,7] [,8]
## [1,] 0.00000000 0.000000000
## [2,] 0.00000000 0.000000000
## [3,] 0.00000000 0.000000000
## [4,] 0.00000000 0.000000000
## [5,] 0.00000000 0.000000000
## [6,] 0.00000000 0.000000000
## [7,] 0.01380088 0.000000000
## [8,] 0.00000000 0.002376322
cov_1f = sigF*beta_%*%t(beta_)+sigeps
cov_1f # This is covariance matrix computed from single factor model (from lm())
## CERN EQIX HFC JKHY KSU
## CERN 3.897737e-03 -8.791295e-06 -3.404636e-06 -7.324709e-07 -3.913437e-06
## EQIX -8.791295e-06 5.703519e-03 7.496920e-05 1.612882e-05 8.617288e-05
## HFC -3.404636e-06 7.496920e-05 1.033759e-02 6.246265e-06 3.337247e-05
## JKHY -7.324709e-07 1.612882e-05 6.246265e-06 1.772552e-03 7.179730e-06
## KSU -3.913437e-06 8.617288e-05 3.337247e-05 7.179730e-06 5.689842e-03
## NTES 1.092253e-06 -2.405113e-05 -9.314366e-06 -2.003886e-06 -1.070634e-05
## NVDA -1.570572e-05 3.458359e-04 1.339331e-04 2.881427e-05 1.539485e-04
## SBAC -6.830641e-06 1.504090e-04 5.824941e-05 1.253173e-05 6.695441e-05
## NTES NVDA SBAC
## CERN 1.092253e-06 -1.570572e-05 -6.830641e-06
## EQIX -2.405113e-05 3.458359e-04 1.504090e-04
## HFC -9.314366e-06 1.339331e-04 5.824941e-05
## JKHY -2.003886e-06 2.881427e-05 1.253173e-05
## KSU -1.070634e-05 1.539485e-04 6.695441e-05
## NTES 8.817481e-03 -4.296754e-05 -1.868719e-05
## NVDA -4.296754e-05 1.441872e-02 2.687068e-04
## SBAC -1.868719e-05 2.687068e-04 2.493187e-03
=======================================================================
Backtesting portfolio using SIT package
=======================================================================
dates <- '2005-01::2021-03'
X <- monthly.ret[dates]
head(X)
## CERN EQIX HFC JKHY KSU
## 2005-01-31 -0.064321874 -0.01871805 0.085755326 0.04419880 -0.015228971
## 2005-02-28 0.046834207 0.04411059 0.245870052 -0.04304979 0.123711670
## 2005-03-31 0.008256439 -0.03311266 -0.009247788 -0.09370291 -0.018348544
## 2005-04-29 0.105503646 -0.17383080 -0.080493556 -0.04446934 -0.017653324
## 2005-05-31 0.125753635 0.08919395 0.115553438 0.03170599 0.056553937
## 2005-06-30 0.040092075 0.13753285 0.223531525 0.03504789 0.009504934
## NTES NVDA SBAC
## 2005-01-31 -0.20351472 -0.02716471 -0.07758616
## 2005-02-28 0.00711766 0.26483441 0.01635514
## 2005-03-31 0.13568941 -0.18040737 0.05057478
## 2005-04-29 0.02447587 -0.07702011 -0.07221048
## 2005-05-31 0.04555602 0.23575029 0.32547203
## 2005-06-30 0.10592582 -0.01402225 0.20106786
dim(X)
## [1] 195 8
tail(X)
## CERN EQIX HFC JKHY KSU
## 2020-10-30 -0.03043305 -0.03800675 -0.06088281 -0.088197332 -0.025935920
## 2020-11-30 0.06777003 -0.04237073 0.28327186 0.085059133 0.056943362
## 2020-12-31 0.05161503 0.02348847 0.10517315 0.009723934 0.098891057
## 2021-01-29 0.02076964 0.03609735 0.10096708 -0.106179561 -0.007152258
## 2021-02-26 -0.13693673 -0.12003306 0.34323077 0.025209088 0.047713000
## 2021-03-31 0.03962966 0.04820003 -0.05543828 0.025298265 0.246255333
## NTES NVDA SBAC
## 2020-10-30 -0.04557146 -0.073648422 -0.088263034
## 2020-11-30 0.04124882 0.069211734 -0.009444331
## 2020-12-31 0.06210449 -0.025567475 -0.017584742
## 2021-01-29 0.20068914 -0.004998005 -0.047708523
## 2021-02-26 -0.04469953 0.055793931 -0.050396423
## 2021-03-31 -0.05945701 -0.026369358 0.090481479
#
prices_monthly <- to.monthly(datas$prices, indexAt = "last", OHLC = FALSE) # indexAt = 'endof'
head(prices_monthly)
## CERN EQIX HFC JKHY KSU NTES NVDA
## 2000-08-31 4.656511 420.3326 0.404515 17.37007 8.275572 0.327738 12.156082
## 2000-09-29 5.681096 230.4527 0.394430 16.99795 7.772331 0.233715 12.538951
## 2000-10-31 7.577343 253.1732 0.398413 21.55360 7.772331 0.150437 9.516680
## 2000-11-30 5.979296 113.6034 0.426302 21.03472 7.548665 0.150437 6.202473
## 2000-12-29 5.658157 113.6034 0.603876 24.36888 9.058400 0.131632 5.017973
## 2001-01-31 5.910480 154.1760 0.579721 17.89667 11.299516 0.096709 7.906240
## SBAC
## 2000-08-31 44.02129
## 2000-09-29 41.37014
## 2000-10-31 49.44688
## 2000-11-30 36.74606
## 2000-12-29 40.50699
## 2001-01-31 44.69948
prices_monthly <- prices_monthly[dates]
head(prices_monthly)
## CERN EQIX HFC JKHY KSU NTES NVDA
## 2005-01-31 6.086342 34.03233 4.220876 16.81811 15.62070 1.811691 7.020281
## 2005-02-28 6.371391 35.53351 5.258663 16.09409 17.55316 1.824586 8.879493
## 2005-03-31 6.423996 34.35690 5.210032 14.58603 17.23109 2.072163 7.277567
## 2005-04-29 7.101751 28.38461 4.790658 13.93740 16.92690 2.122881 6.717048
## 2005-05-31 7.994822 30.91635 5.344235 14.37930 17.88419 2.219591 8.300594
## 2005-06-30 8.315351 35.16837 6.538840 14.88326 18.05417 2.454703 8.184201
## SBAC
## 2005-01-31 8.444196
## 2005-02-28 8.582302
## 2005-03-31 9.016350
## 2005-04-29 8.365275
## 2005-05-31 11.087938
## 2005-06-30 13.317366
tail(prices_monthly)
## CERN EQIX HFC JKHY KSU NTES NVDA
## 2020-10-30 69.68066 722.6181 18.06076 146.9587 174.9587 86.32986 500.9242
## 2020-11-30 74.40292 692.0002 23.17686 159.4589 184.9215 89.89086 535.5940
## 2020-12-31 78.24322 708.2543 25.61445 161.0095 203.2085 95.47349 521.9002
## 2021-01-29 79.86831 733.8204 28.20067 143.9136 201.7551 114.63398 519.2917
## 2021-02-26 68.93140 645.7377 37.88000 147.5415 211.3815 109.50990 548.2651
## 2021-03-31 71.66313 676.8622 35.78000 151.2740 263.4353 102.99876 533.8077
## SBAC
## 2020-10-30 288.6498
## 2020-11-30 285.9237
## 2020-12-31 280.8958
## 2021-01-29 267.4947
## 2021-02-26 254.0139
## 2021-03-31 276.9974
#*****************************************************************
# Code Strategies
#******************************************************************
# We need to create environment variable required by SIT
# data_m <- new.env()
#
# data_m$KSU <- prices_monthly$KSU
# colnames(data_m$KSU) <- 'Close'
# head(data_m$KSU)
# Change column names to Close which will be consistent with the requirement of SIT package
# data_m$HFC <- prices_monthly$HFC %>% `colnames<-` (c("Close"))
# data_m$KSU <- prices_monthly$KSU %>% `colnames<-` (c("Close"))
# data_m$JKHY <- prices_monthly$JKHY %>% `colnames<-` (c("Close"))
# data_m$CERN <- prices_monthly$CERN %>% `colnames<-` (c("Close"))
# data_m$NTES <- prices_monthly$NTES %>% `colnames<-` (c("Close"))
# data_m$NVDA <- prices_monthly$NVDA %>% `colnames<-` (c("Close"))
# data_m$SBAC <- prices_monthly$SBAC %>% `colnames<-` (c("Close"))
# data_m$EQIX <- prices_monthly$EQIX %>% `colnames<-` (c("Close"))
# names(data_m)
# head(data_m$EQIX)
Using loop to save time in naming
data_m <- list()
#i = 1
for ( i in 1:length(tickers)){
data_m[[tickers[i]]] <- prices_monthly[,i] %>% `colnames<-` (c("Close"))
}
# data_env <- new.env()
# list2env(data_m, envir = data_env)
# convert from list to environment variable because SIT package requires
# the input data to be environmental variable which can be processed using bt.prep()
data_m <- list2env(data_m)
names(data_m)
## [1] "EQIX" "SBAC" "NVDA" "NTES" "CERN" "JKHY" "KSU" "HFC"
Check if the column name is ‘Close’
head(data_m$NVDA)
## Close
## 2005-01-31 1.811691
## 2005-02-28 1.824586
## 2005-03-31 2.072163
## 2005-04-29 2.122881
## 2005-05-31 2.219591
## 2005-06-30 2.454703
bt.prep(data_m, align='remove.na', dates = dates)
names(data_m)
## [1] "prices" "execution.price" "weight" "dates"
## [5] "symbolnames" "EQIX" "SBAC" "NVDA"
## [9] "NTES" "CERN" "JKHY" "KSU"
## [13] "HFC"
===============================================================
Equal weighting portfolio
Equal Weight 1/N Benchmark
===============================================================
models <- list()
prices <- data_m$prices
# data_m$weight[] = NA
N <- length(tickers)
data_m$weight = ntop(prices, N)
head(data_m$weight)
## CERN EQIX HFC JKHY KSU NTES NVDA SBAC
## 2005-01-31 08:00:00 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## 2005-02-28 08:00:00 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## 2005-03-31 08:00:00 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## 2005-04-29 08:00:00 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## 2005-05-31 08:00:00 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## 2005-06-30 08:00:00 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
data_m$weight[1:59, ] <- NA
models$equal.weight = bt.run.share(data_m)
## Latest weights :
## CERN EQIX HFC JKHY KSU NTES NVDA SBAC
## 2021-03-31 08:00:00 8.54 8.69 3.94 5.7 8.64 8.18 19.51 36.81
##
## Performance summary :
## CAGR Best Worst
## 16 13.2 -12.1
# head(models$equal.weight$ret, 62)
# Slightly difference between bt.run.share() and bt.run()
capital = 100000
data_m$weight[] = (capital / prices) * data_m$weight
models$equal.weight.share = bt.run(data_m, type='share')
## Latest weights :
## CERN EQIX HFC JKHY KSU NTES NVDA SBAC
## 2021-03-31 08:00:00 12.5 12.5 12.5 12.5 12.5 12.5 12.5 12.5
##
## Performance summary :
## CAGR Best Worst
## 16.9 13.6 -9.6
# head(models$equal.weight.share$ret, 62)
# head(equal.weight$ret)
================================================================
MVP portfolio
*****************************************************************
Create Constraints
*****************************************************************
constraints = new.constraints(N, lb = -Inf, ub = +Inf)
# SUM x.i = 1
constraints = add.constraints(rep(1, N), 1, type = '=', constraints)
#
ret = prices / mlag(prices) - 1
weight = coredata(ret)
weight[] = NA
#
# head(data_m$weight)
# compute covariance matrix based on historical 60 months returns
# i = 60
# To make covariance matrix estimate more stable, use the Ledoit-Wolf covariance shrinkage estimator from tawny package
# 1. ia$cov = tawny::cov.shrink(hist) or
# 2. ia$cov = cor(coredata(hist), use='complete.obs', method='spearman') * (s0 %*% t(s0)) or
# 3. ia$cov = cor(coredata(hist), use='complete.obs', method='kendall') * (s0 %*% t(s0)
# i = 60
for( i in 60:dim(weight)[1]) {
hist = ret[ (i- 60 + 1):i, ]
# create historical input assumptions
ia = create.historical.ia(hist, 12) # 12 is annulized factor for monthly returns
s0 = apply(na.omit(coredata(hist)), 2, sd)
ia$cov = cor(coredata(hist), use='complete.obs',method='pearson') * (s0 %*% t(s0))
weight[i,] = min.risk.portfolio(ia, constraints) # use SIT's function min.risk.portfolio()
}
## Loading required package: kernlab
##
## Attaching package: 'kernlab'
## The following object is masked from 'package:SIT':
##
## cross
## The following object is masked from 'package:purrr':
##
## cross
dim(weight)
## [1] 195 8
# 195 8
head(weight, 70)
## CERN EQIX HFC JKHY KSU NTES
## [1,] NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA
## [24,] NA NA NA NA NA NA
## [25,] NA NA NA NA NA NA
## [26,] NA NA NA NA NA NA
## [27,] NA NA NA NA NA NA
## [28,] NA NA NA NA NA NA
## [29,] NA NA NA NA NA NA
## [30,] NA NA NA NA NA NA
## [31,] NA NA NA NA NA NA
## [32,] NA NA NA NA NA NA
## [33,] NA NA NA NA NA NA
## [34,] NA NA NA NA NA NA
## [35,] NA NA NA NA NA NA
## [36,] NA NA NA NA NA NA
## [37,] NA NA NA NA NA NA
## [38,] NA NA NA NA NA NA
## [39,] NA NA NA NA NA NA
## [40,] NA NA NA NA NA NA
## [41,] NA NA NA NA NA NA
## [42,] NA NA NA NA NA NA
## [43,] NA NA NA NA NA NA
## [44,] NA NA NA NA NA NA
## [45,] NA NA NA NA NA NA
## [46,] NA NA NA NA NA NA
## [47,] NA NA NA NA NA NA
## [48,] NA NA NA NA NA NA
## [49,] NA NA NA NA NA NA
## [50,] NA NA NA NA NA NA
## [51,] NA NA NA NA NA NA
## [52,] NA NA NA NA NA NA
## [53,] NA NA NA NA NA NA
## [54,] NA NA NA NA NA NA
## [55,] NA NA NA NA NA NA
## [56,] NA NA NA NA NA NA
## [57,] NA NA NA NA NA NA
## [58,] NA NA NA NA NA NA
## [59,] NA NA NA NA NA NA
## [60,] 0.6506858 0.1630486 0.1640682 0.09980743 -0.2365001 0.037552221
## [61,] 0.6613121 0.1667538 0.1620214 0.10792216 -0.2346133 0.030297899
## [62,] 0.6721412 0.1666004 0.1571951 0.11463660 -0.2156373 0.017145468
## [63,] 0.6790111 0.1628983 0.1525422 0.11138608 -0.2152667 0.016993072
## [64,] 0.6783518 0.1645960 0.1490279 0.11242843 -0.2124555 0.014688017
## [65,] 0.6656822 0.2056694 0.1621668 0.10819839 -0.2300020 0.003052561
## [66,] 0.6582386 0.2175427 0.1573987 0.11866650 -0.2162170 -0.008461939
## [67,] 0.6500353 0.2563896 0.1712984 0.10415797 -0.2460031 -0.005046055
## [68,] 0.6488025 0.2440374 0.1709201 0.11895234 -0.2785601 -0.015995144
## [69,] 0.6397226 0.2345371 0.1693333 0.12301158 -0.2830282 -0.009728167
## [70,] 0.6324751 0.2436580 0.1559293 0.11847589 -0.2135562 -0.032412507
## NVDA SBAC
## [1,] NA NA
## [2,] NA NA
## [3,] NA NA
## [4,] NA NA
## [5,] NA NA
## [6,] NA NA
## [7,] NA NA
## [8,] NA NA
## [9,] NA NA
## [10,] NA NA
## [11,] NA NA
## [12,] NA NA
## [13,] NA NA
## [14,] NA NA
## [15,] NA NA
## [16,] NA NA
## [17,] NA NA
## [18,] NA NA
## [19,] NA NA
## [20,] NA NA
## [21,] NA NA
## [22,] NA NA
## [23,] NA NA
## [24,] NA NA
## [25,] NA NA
## [26,] NA NA
## [27,] NA NA
## [28,] NA NA
## [29,] NA NA
## [30,] NA NA
## [31,] NA NA
## [32,] NA NA
## [33,] NA NA
## [34,] NA NA
## [35,] NA NA
## [36,] NA NA
## [37,] NA NA
## [38,] NA NA
## [39,] NA NA
## [40,] NA NA
## [41,] NA NA
## [42,] NA NA
## [43,] NA NA
## [44,] NA NA
## [45,] NA NA
## [46,] NA NA
## [47,] NA NA
## [48,] NA NA
## [49,] NA NA
## [50,] NA NA
## [51,] NA NA
## [52,] NA NA
## [53,] NA NA
## [54,] NA NA
## [55,] NA NA
## [56,] NA NA
## [57,] NA NA
## [58,] NA NA
## [59,] NA NA
## [60,] 0.11143070 0.009907210
## [61,] 0.10453653 0.001769320
## [62,] 0.09188497 -0.003966504
## [63,] 0.08831608 0.004119795
## [64,] 0.08773583 0.005627491
## [65,] 0.07656245 0.008670281
## [66,] 0.07410197 -0.001269515
## [67,] 0.06785439 0.001313592
## [68,] 0.09515242 0.016690554
## [69,] 0.11792703 0.008224775
## [70,] 0.10556000 -0.010129462
tail(weight)
## CERN EQIX HFC JKHY KSU NTES NVDA
## [190,] 0.3658217 0.1230770 0.1331788 0.04368771 0.2283614 0.09483016 0.05911254
## [191,] 0.3475902 0.1438089 0.1345830 0.03820816 0.2402433 0.08958937 0.05029831
## [192,] 0.3457549 0.1535327 0.1199478 0.03758611 0.2318964 0.11210127 0.04394611
## [193,] 0.3060546 0.1563895 0.1221102 0.04764332 0.2346858 0.10839713 0.06097096
## [194,] 0.3068238 0.1731704 0.1238168 0.06548481 0.2110699 0.07688665 0.07532633
## [195,] 0.2992969 0.1690535 0.1336878 0.07556138 0.2252926 0.03486647 0.08457208
## SBAC
## [190,] -0.04806936
## [191,] -0.04432127
## [192,] -0.04476527
## [193,] -0.03625162
## [194,] -0.03257867
## [195,] -0.02233076
Assign minimum variance weights to data_m$weight
capital = 100000
data_m$weight <- data_m$prices
data_m$weight[] <- NA
data_m$weight[] <- weight
data_m$weight[] = (capital / prices) * data_m$weight
models$mvp.hist.cov = bt.run(data_m, type='share')
## Latest weights :
## CERN EQIX HFC JKHY KSU NTES NVDA SBAC
## 2021-03-31 08:00:00 30.68 17.32 12.38 6.55 21.11 7.69 7.53 -3.26
##
## Performance summary :
## CAGR Best Worst
## 14.9 12.3 -9.4
=================================================================
CAPM model (single factor model)
*****************************************************************
Create Constraints
*****************************************************************
constraints = new.constraints(N, lb = -Inf, ub = +Inf)
# SUM x.i = 1
constraints = add.constraints(rep(1, N), 1, type = '=', constraints)
ret = prices / mlag(prices) - 1
weight_1 = coredata(ret)
weight_1[] = NA
#
# head(data_m$weight)
# compute covariance matrix based on historical 60 months returns
# i = 60
f <- fama_lib_month[dates]/100
for( i in 60:dim(weight_1)[1]) {
hist = ret[ (i- 60 + 1):i, ]
Xi <- hist
fi <- f$Mkt.RF[(i - 60 + 1):i, ]
fiti = lm(formula = Xi ~ fi)
sigF = as.numeric(var(fi))
beta_ = as.matrix(fiti$coefficients)
beta_ = as.matrix(beta_[-1,])
sigeps = crossprod(fiti$residuals)/(T-2)
# sigeps = as.matrix(var(fit$residuals)) # you can use this way too
sigeps = diag(diag(sigeps))
cov_1f = sigF*beta_%*%t(beta_)+sigeps
# cov_1f
ia$cov = cov_1f
weight_1[i, ] = min.risk.portfolio(ia, constraints)
}
dim(weight_1)
## [1] 195 8
# 195 8
head(weight_1, 70)
## CERN EQIX HFC JKHY KSU NTES
## [1,] NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA
## [24,] NA NA NA NA NA NA
## [25,] NA NA NA NA NA NA
## [26,] NA NA NA NA NA NA
## [27,] NA NA NA NA NA NA
## [28,] NA NA NA NA NA NA
## [29,] NA NA NA NA NA NA
## [30,] NA NA NA NA NA NA
## [31,] NA NA NA NA NA NA
## [32,] NA NA NA NA NA NA
## [33,] NA NA NA NA NA NA
## [34,] NA NA NA NA NA NA
## [35,] NA NA NA NA NA NA
## [36,] NA NA NA NA NA NA
## [37,] NA NA NA NA NA NA
## [38,] NA NA NA NA NA NA
## [39,] NA NA NA NA NA NA
## [40,] NA NA NA NA NA NA
## [41,] NA NA NA NA NA NA
## [42,] NA NA NA NA NA NA
## [43,] NA NA NA NA NA NA
## [44,] NA NA NA NA NA NA
## [45,] NA NA NA NA NA NA
## [46,] NA NA NA NA NA NA
## [47,] NA NA NA NA NA NA
## [48,] NA NA NA NA NA NA
## [49,] NA NA NA NA NA NA
## [50,] NA NA NA NA NA NA
## [51,] NA NA NA NA NA NA
## [52,] NA NA NA NA NA NA
## [53,] NA NA NA NA NA NA
## [54,] NA NA NA NA NA NA
## [55,] NA NA NA NA NA NA
## [56,] NA NA NA NA NA NA
## [57,] NA NA NA NA NA NA
## [58,] NA NA NA NA NA NA
## [59,] NA NA NA NA NA NA
## [60,] 0.3515628 0.09776704 0.1686894 0.06928919 0.08905936 0.08267310
## [61,] 0.3478773 0.10174204 0.1661870 0.07141374 0.08914277 0.08360828
## [62,] 0.3470957 0.10141917 0.1642139 0.07555324 0.08869669 0.08287205
## [63,] 0.3522214 0.10032255 0.1627082 0.07425723 0.08807766 0.08196803
## [64,] 0.3465917 0.10205359 0.1619331 0.07550750 0.09252555 0.08086526
## [65,] 0.3352742 0.11730658 0.1618044 0.07615146 0.09075831 0.08044314
## [66,] 0.3317393 0.12498348 0.1570592 0.08021834 0.09021584 0.08025179
## [67,] 0.3278800 0.13401093 0.1586723 0.07973053 0.08770203 0.08074009
## [68,] 0.3205098 0.13470513 0.1573716 0.08227631 0.09049018 0.08086428
## [69,] 0.3178098 0.13229018 0.1560407 0.08329108 0.09087460 0.08239456
## [70,] 0.3187709 0.13519500 0.1537719 0.08392917 0.08791237 0.08125975
## NVDA SBAC
## [1,] NA NA
## [2,] NA NA
## [3,] NA NA
## [4,] NA NA
## [5,] NA NA
## [6,] NA NA
## [7,] NA NA
## [8,] NA NA
## [9,] NA NA
## [10,] NA NA
## [11,] NA NA
## [12,] NA NA
## [13,] NA NA
## [14,] NA NA
## [15,] NA NA
## [16,] NA NA
## [17,] NA NA
## [18,] NA NA
## [19,] NA NA
## [20,] NA NA
## [21,] NA NA
## [22,] NA NA
## [23,] NA NA
## [24,] NA NA
## [25,] NA NA
## [26,] NA NA
## [27,] NA NA
## [28,] NA NA
## [29,] NA NA
## [30,] NA NA
## [31,] NA NA
## [32,] NA NA
## [33,] NA NA
## [34,] NA NA
## [35,] NA NA
## [36,] NA NA
## [37,] NA NA
## [38,] NA NA
## [39,] NA NA
## [40,] NA NA
## [41,] NA NA
## [42,] NA NA
## [43,] NA NA
## [44,] NA NA
## [45,] NA NA
## [46,] NA NA
## [47,] NA NA
## [48,] NA NA
## [49,] NA NA
## [50,] NA NA
## [51,] NA NA
## [52,] NA NA
## [53,] NA NA
## [54,] NA NA
## [55,] NA NA
## [56,] NA NA
## [57,] NA NA
## [58,] NA NA
## [59,] NA NA
## [60,] 0.08862724 0.05233182
## [61,] 0.08745864 0.05257020
## [62,] 0.08501849 0.05513075
## [63,] 0.08431607 0.05612892
## [64,] 0.08367287 0.05685041
## [65,] 0.08069680 0.05756503
## [66,] 0.08014636 0.05538567
## [67,] 0.07689143 0.05437263
## [68,] 0.07930542 0.05447727
## [69,] 0.08420018 0.05309897
## [70,] 0.08615189 0.05300901
tail(weight_1)
## CERN EQIX HFC JKHY KSU NTES NVDA
## [190,] 0.2214719 0.1635502 0.1618538 0.03390453 0.2049068 0.1173410 0.05643956
## [191,] 0.2125868 0.1762334 0.1662881 0.03170385 0.1979474 0.1182107 0.05732418
## [192,] 0.2110421 0.1733044 0.1612622 0.03211930 0.1949625 0.1317271 0.05671826
## [193,] 0.1978185 0.1717154 0.1703870 0.03333870 0.2030716 0.1311372 0.05369844
## [194,] 0.1998078 0.1733260 0.1706336 0.03150867 0.1869612 0.1410736 0.05753326
## [195,] 0.2054456 0.1737738 0.1747896 0.03272071 0.1968290 0.1165173 0.05914727
## SBAC
## [190,] 0.04053226
## [191,] 0.03970555
## [192,] 0.03886415
## [193,] 0.03883332
## [194,] 0.03915590
## [195,] 0.04077671
Assign minimum variance weights to data_m$weight
capital = 100000
data_m$weight <- data_m$prices
data_m$weight[] <- NA
data_m$weight[] <- weight_1
data_m$weight[] = (capital / prices) * data_m$weight
models$mvp.capm.cov = bt.run(data_m, type='share')
## Latest weights :
## CERN EQIX HFC JKHY KSU NTES NVDA SBAC
## 2021-03-31 08:00:00 19.98 17.33 17.06 3.15 18.7 14.11 5.75 3.92
##
## Performance summary :
## CAGR Best Worst
## 15.5 11.3 -8.6
models$mvp.capm.cov$cagr
## [1] 0.1552803
models$mvp.hist.cov$cagr
## [1] 0.1488864
=================================================================
FF 3 factor model
*****************************************************************
Create Constraints
*****************************************************************
constraints = new.constraints(N, lb = -Inf, ub = +Inf)
# SUM x.i = 1
constraints = add.constraints(rep(1, N), 1, type = '=', constraints)
ret = prices / mlag(prices) - 1
weight_3 = coredata(ret)
weight_3[] = NA
#
# head(data_m$weight)
# compute covariance matrix based on historical 60 months returns
# i = 60
f <- fama_lib_month[dates]/100
for( i in 60:dim(weight_1)[1]) {
hist = ret[ (i- 60 + 1):i, ]
Xi <- hist
fi <- f[(i - 60 + 1):i, ]
fiti = lm(formula = Xi ~ fi)
sigF = as.matrix(var(fi))
beta_ = as.matrix(fiti$coefficients)
beta_ = as.matrix(beta_[-1,])
sigeps = crossprod(fiti$residuals)/(T-4) # note (T - 4)
# sigeps = as.matrix(var(fit$residuals)) # you can use this way too
sigeps = diag(diag(sigeps))
cov_3f = t(beta_)%*% sigF %*% beta_ + sigeps
# cov_1f
ia$cov = cov_3f
weight_3[i, ] = min.risk.portfolio(ia, constraints)
}
dim(weight_3)
## [1] 195 8
# 195 8
head(weight_3, 70)
## CERN EQIX HFC JKHY KSU NTES NVDA
## [1,] NA NA NA NA NA NA NA
## [2,] NA NA NA NA NA NA NA
## [3,] NA NA NA NA NA NA NA
## [4,] NA NA NA NA NA NA NA
## [5,] NA NA NA NA NA NA NA
## [6,] NA NA NA NA NA NA NA
## [7,] NA NA NA NA NA NA NA
## [8,] NA NA NA NA NA NA NA
## [9,] NA NA NA NA NA NA NA
## [10,] NA NA NA NA NA NA NA
## [11,] NA NA NA NA NA NA NA
## [12,] NA NA NA NA NA NA NA
## [13,] NA NA NA NA NA NA NA
## [14,] NA NA NA NA NA NA NA
## [15,] NA NA NA NA NA NA NA
## [16,] NA NA NA NA NA NA NA
## [17,] NA NA NA NA NA NA NA
## [18,] NA NA NA NA NA NA NA
## [19,] NA NA NA NA NA NA NA
## [20,] NA NA NA NA NA NA NA
## [21,] NA NA NA NA NA NA NA
## [22,] NA NA NA NA NA NA NA
## [23,] NA NA NA NA NA NA NA
## [24,] NA NA NA NA NA NA NA
## [25,] NA NA NA NA NA NA NA
## [26,] NA NA NA NA NA NA NA
## [27,] NA NA NA NA NA NA NA
## [28,] NA NA NA NA NA NA NA
## [29,] NA NA NA NA NA NA NA
## [30,] NA NA NA NA NA NA NA
## [31,] NA NA NA NA NA NA NA
## [32,] NA NA NA NA NA NA NA
## [33,] NA NA NA NA NA NA NA
## [34,] NA NA NA NA NA NA NA
## [35,] NA NA NA NA NA NA NA
## [36,] NA NA NA NA NA NA NA
## [37,] NA NA NA NA NA NA NA
## [38,] NA NA NA NA NA NA NA
## [39,] NA NA NA NA NA NA NA
## [40,] NA NA NA NA NA NA NA
## [41,] NA NA NA NA NA NA NA
## [42,] NA NA NA NA NA NA NA
## [43,] NA NA NA NA NA NA NA
## [44,] NA NA NA NA NA NA NA
## [45,] NA NA NA NA NA NA NA
## [46,] NA NA NA NA NA NA NA
## [47,] NA NA NA NA NA NA NA
## [48,] NA NA NA NA NA NA NA
## [49,] NA NA NA NA NA NA NA
## [50,] NA NA NA NA NA NA NA
## [51,] NA NA NA NA NA NA NA
## [52,] NA NA NA NA NA NA NA
## [53,] NA NA NA NA NA NA NA
## [54,] NA NA NA NA NA NA NA
## [55,] NA NA NA NA NA NA NA
## [56,] NA NA NA NA NA NA NA
## [57,] NA NA NA NA NA NA NA
## [58,] NA NA NA NA NA NA NA
## [59,] NA NA NA NA NA NA NA
## [60,] 0.3607411 0.1047147 0.1569195 0.05453581 0.08995028 0.07068614 0.10481149
## [61,] 0.3566157 0.1074269 0.1554380 0.05707153 0.09019968 0.07280581 0.10263258
## [62,] 0.3561447 0.1082406 0.1534349 0.06012103 0.08944039 0.07201198 0.10077301
## [63,] 0.3620940 0.1062519 0.1518872 0.05971436 0.08873464 0.07145021 0.09780437
## [64,] 0.3579641 0.1093819 0.1510569 0.05955031 0.09063873 0.06928813 0.09887286
## [65,] 0.3475084 0.1234296 0.1498040 0.05966458 0.08918756 0.06901911 0.09762726
## [66,] 0.3466258 0.1305142 0.1446279 0.06089690 0.08808157 0.06898158 0.09779115
## [67,] 0.3449019 0.1389263 0.1452142 0.05960664 0.08500881 0.06850396 0.09413168
## [68,] 0.3361219 0.1379006 0.1431892 0.06418375 0.08476967 0.06517611 0.10201902
## [69,] 0.3272498 0.1341506 0.1437109 0.06775499 0.08536517 0.06905234 0.10772096
## [70,] 0.3309682 0.1361130 0.1413399 0.06559347 0.08318418 0.06683858 0.10900390
## SBAC
## [1,] NA
## [2,] NA
## [3,] NA
## [4,] NA
## [5,] NA
## [6,] NA
## [7,] NA
## [8,] NA
## [9,] NA
## [10,] NA
## [11,] NA
## [12,] NA
## [13,] NA
## [14,] NA
## [15,] NA
## [16,] NA
## [17,] NA
## [18,] NA
## [19,] NA
## [20,] NA
## [21,] NA
## [22,] NA
## [23,] NA
## [24,] NA
## [25,] NA
## [26,] NA
## [27,] NA
## [28,] NA
## [29,] NA
## [30,] NA
## [31,] NA
## [32,] NA
## [33,] NA
## [34,] NA
## [35,] NA
## [36,] NA
## [37,] NA
## [38,] NA
## [39,] NA
## [40,] NA
## [41,] NA
## [42,] NA
## [43,] NA
## [44,] NA
## [45,] NA
## [46,] NA
## [47,] NA
## [48,] NA
## [49,] NA
## [50,] NA
## [51,] NA
## [52,] NA
## [53,] NA
## [54,] NA
## [55,] NA
## [56,] NA
## [57,] NA
## [58,] NA
## [59,] NA
## [60,] 0.05764096
## [61,] 0.05780987
## [62,] 0.05983344
## [63,] 0.06206337
## [64,] 0.06324704
## [65,] 0.06375944
## [66,] 0.06248095
## [67,] 0.06370650
## [68,] 0.06663974
## [69,] 0.06499519
## [70,] 0.06695875
tail(weight_3)
## CERN EQIX HFC JKHY KSU NTES NVDA
## [190,] 0.2535843 0.1420366 0.1384924 0.04276513 0.2125199 0.11899815 0.04770770
## [191,] 0.2329274 0.1654764 0.1440556 0.03308800 0.2203344 0.11612072 0.04542731
## [192,] 0.2306288 0.1596320 0.1384554 0.03379871 0.2172846 0.13502086 0.04352262
## [193,] 0.2040688 0.1587935 0.1494872 0.03699029 0.2232434 0.13657466 0.04840729
## [194,] 0.2027429 0.1646556 0.1584125 0.02741430 0.2101449 0.13823556 0.05532834
## [195,] 0.2109423 0.1653614 0.1648384 0.03113387 0.2201402 0.09856064 0.06242806
## SBAC
## [190,] 0.04389580
## [191,] 0.04257018
## [192,] 0.04165699
## [193,] 0.04243491
## [194,] 0.04306587
## [195,] 0.04659506
Assign minimum variance weights to data_m$weight
capital = 100000
data_m$weight <- data_m$prices
data_m$weight[] <- NA
data_m$weight[] <- weight_3
data_m$weight[] = (capital / prices) * data_m$weight
models$mvp.ff3.cov = bt.run(data_m, type='share')
## Latest weights :
## CERN EQIX HFC JKHY KSU NTES NVDA SBAC
## 2021-03-31 08:00:00 20.27 16.47 15.84 2.74 21.01 13.82 5.53 4.31
##
## Performance summary :
## CAGR Best Worst
## 15.6 11.4 -9.1
=============================================================
Principal Component Analysis (PCA)
=============================================================
library(covFactorModel)
X_i <-X[1:60,]
factor_pca <- factorModel(X_i, type = "S", K = 3, max_iter = 10)
cbind(alpha = factor_pca$alpha, beta = factor_pca$beta)
## alpha factor1 factor2 factor3
## CERN 0.023270140 -0.04987022 0.005325561 -0.008727843
## EQIX 0.022552077 -0.09904223 0.030680877 -0.016892567
## HFC 0.020222607 -0.05493627 0.007170780 0.023949692
## JKHY 0.005672384 -0.02687649 0.011404884 -0.007012624
## KSU 0.018269879 -0.05319466 0.039958217 0.040898599
## NTES 0.025794575 -0.06133519 0.034139829 -0.053821596
## NVDA 0.026442881 -0.12675785 -0.081364792 0.001852985
## SBAC 0.028204138 -0.06876748 0.030382689 0.027219197
#
# Statistical 3-factor model
K <- 3
X_trn <- X_i
T_trn <- dim(X_i)[1]
alpha <- colMeans(X_trn)
X_trn_ <- X_trn - matrix(alpha, T_trn, N, byrow = TRUE)
Sigma_prev <- matrix(0, N, N)
Sigma <- (1/(T_trn-1)) * t(X_trn_) %*% X_trn_
eigSigma <- eigen(Sigma)
while (norm(Sigma - Sigma_prev, "F")/norm(Sigma, "F") > 1e-3) {
B <- eigSigma$vectors[, 1:K] %*% diag(sqrt(eigSigma$values[1:K]), K, K)
Psi <- diag(diag(Sigma - B %*% t(B)))
Sigma_prev <- Sigma
Sigma <- B %*% t(B) + Psi
eigSigma <- eigen(Sigma - Psi)
}
B: factor loadings
B
## [,1] [,2] [,3]
## [1,] -0.05361747 0.018608974 -0.002152947
## [2,] -0.10268075 0.020996577 -0.024907970
## [3,] -0.07336302 -0.075841600 0.027041843
## [4,] -0.02717521 0.008802363 -0.011196077
## [5,] -0.06354532 -0.063424136 -0.041483092
## [6,] -0.07199784 0.052932895 -0.069369422
## [7,] -0.11830733 0.038915744 0.083830900
## [8,] -0.07637508 -0.028984746 -0.016941602
# fiti = lm(formula = X_i ~ t(B))
Sigma_PCA3 <- Sigma
Sigma_PCA3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.008804419 0.005949832 0.0024639851 0.001644973 0.002316192
## [2,] 0.005949832 0.015092330 0.0052669985 0.003254062 0.006226451
## [3,] 0.002463985 0.005266999 0.0175598543 0.001023307 0.008350286
## [4,] 0.001644973 0.003254062 0.0010233073 0.004328908 0.001633023
## [5,] 0.002316192 0.006226451 0.0083502855 0.001633023 0.015041805
## [6,] 0.004994717 0.010232053 -0.0006084137 0.003199156 0.004095561
## [7,] 0.006887038 0.010876926 0.0079948929 0.002618999 0.001572115
## [8,] 0.003592136 0.007655651 0.0073432235 0.002010054 0.007394401
## [,6] [,7] [,8]
## [1,] 0.0049947171 0.006887038 0.003592136
## [2,] 0.0102320532 0.010876926 0.007655651
## [3,] -0.0006084137 0.007994893 0.007343223
## [4,] 0.0031991560 0.002618999 0.002010054
## [5,] 0.0040955609 0.001572115 0.007394401
## [6,] 0.0170421169 0.004762494 0.005139823
## [7,] 0.0047624940 0.023632553 0.006487539
## [8,] 0.0051398229 0.006487539 0.012804919
diag(Sigma_PCA3)
## [1] 0.008804419 0.015092330 0.017559854 0.004328908 0.015041805 0.017042117
## [7] 0.023632553 0.012804919
error
norm(Sigma_PCA3 - cov(X_i), "F")
## [1] 0.01146673
—————————————————————————
By Eric Zivot
use R princomp() function for principal component analysis
—————————————————————————
pc.fit = princomp(X_i)
class(pc.fit)
## [1] "princomp"
names(pc.fit)
## [1] "sdev" "loadings" "center" "scale" "n.obs" "scores" "call"
pc.fit
## Call:
## princomp(x = X_i)
##
## Standard deviations:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## 0.21858735 0.12462898 0.12277528 0.11291491 0.08753967 0.07798872 0.07096955
## Comp.8
## 0.04983295
##
## 8 variables and 60 observations.
summary(pc.fit)
## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 0.2185874 0.1246290 0.1227753 0.1129149 0.08753967
## Proportion of Variance 0.4250860 0.1381863 0.1341061 0.1134304 0.06817679
## Cumulative Proportion 0.4250860 0.5632723 0.6973784 0.8108088 0.87898556
## Comp.6 Comp.7 Comp.8
## Standard deviation 0.07798872 0.07096955 0.04983295
## Proportion of Variance 0.05411159 0.04480958 0.02209327
## Cumulative Proportion 0.93309715 0.97790673 1.00000000
#eigenvalues & eigenvectors
eig.value <- eigen(cov(X_i))$values
pc.fit
## Call:
## princomp(x = X_i)
##
## Standard deviations:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## 0.21858735 0.12462898 0.12277528 0.11291491 0.08753967 0.07798872 0.07096955
## Comp.8
## 0.04983295
##
## 8 variables and 60 observations.
eig.vec <- eigen(cov(X_i))$vectors
#
B_hat <- eig.vec[, 1:3]%*%diag(sqrt(eig.value[1:3]), 3, 3)
B_hat
## [,1] [,2] [,3]
## [1,] -0.05361747 0.018608974 0.002152947
## [2,] -0.10268075 0.020996577 0.024907970
## [3,] -0.07336302 -0.075841600 -0.027041843
## [4,] -0.02717521 0.008802363 0.011196077
## [5,] -0.06354532 -0.063424136 0.041483092
## [6,] -0.07199784 0.052932895 0.069369422
## [7,] -0.11830733 0.038915744 -0.083830900
## [8,] -0.07637508 -0.028984746 0.016941602
#
#
plot(pc.fit)
loadings(pc.fit)
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## CERN 0.243 0.148 0.153 0.821 0.458 0.102
## EQIX 0.466 0.167 0.201 0.170 0.165 0.117 -0.677 0.433
## HFC 0.333 -0.603 -0.218 -0.637 0.214 -0.138
## JKHY 0.123 0.164 -0.134 0.234 -0.370 -0.857
## KSU 0.288 -0.505 0.335 0.472 -0.502 -0.191 0.184
## NTES 0.327 0.421 0.560 -0.512 -0.193 -0.204 0.231
## NVDA 0.537 0.310 -0.677 0.117 -0.191 -0.290 0.141
## SBAC 0.346 -0.231 0.137 0.127 0.789 -0.221 0.263 -0.228
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8
## SS loadings 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125
## Cumulative Var 0.125 0.250 0.375 0.500 0.625 0.750 0.875 1.000
pc.fit$loadings[, 1:3]
## Comp.1 Comp.2 Comp.3
## CERN 0.2432381 0.14806546 0.01738893
## EQIX 0.4658160 0.16706284 0.20117675
## HFC 0.3328147 -0.60344658 -0.21841162
## JKHY 0.1232816 0.07003749 0.09042850
## KSU 0.2882763 -0.50464492 0.33505073
## NTES 0.3266215 0.42116958 0.56028310
## NVDA 0.5367067 0.30963973 -0.67708560
## SBAC 0.3464791 -0.23062206 0.13683397
# pc factors are in the scores component. Note these scores are based on
# centered data
head(pc.fit$scores[, 1:3])
## Comp.1 Comp.2 Comp.3
## 2005-01-31 -0.1639501 -0.128585797 -0.13853565
## 2005-02-28 0.2330032 -0.117016589 -0.18711190
## 2005-03-31 -0.1295699 -0.005143263 0.17841038
## 2005-04-29 -0.2122847 0.045327324 0.02292385
## 2005-05-31 0.3237324 -0.044139668 -0.08042189
## 2005-06-30 0.1907589 -0.113152694 0.07468680
pc.fit$scores[, 1:3]
## Comp.1 Comp.2 Comp.3
## 2005-01-31 -0.163950105 -0.128585797 -0.138535653
## 2005-02-28 0.233003218 -0.117016589 -0.187111902
## 2005-03-31 -0.129569915 -0.005143263 0.178410377
## 2005-04-29 -0.212284705 0.045327324 0.022923850
## 2005-05-31 0.323732402 -0.044139668 -0.080421891
## 2005-06-30 0.190758905 -0.113152694 0.074686796
## 2005-07-29 0.113915948 -0.072520149 0.082203899
## 2005-08-31 0.047590584 0.082597684 -0.082139640
## 2005-09-30 0.236588793 -0.004566832 0.088521197
## 2005-10-31 -0.253968414 -0.003127236 -0.109384364
## 2005-11-30 0.119930888 -0.196920116 -0.112419644
## 2005-12-30 -0.090741775 0.029640220 -0.024910947
## 2006-01-31 0.405708538 -0.018329632 0.022529139
## 2006-02-28 0.000738157 0.287258355 0.124874542
## 2006-03-31 0.351376242 -0.009826165 -0.059486066
## 2006-04-28 -0.083983704 -0.096198587 -0.090492038
## 2006-05-31 -0.255789725 -0.200447713 0.067471928
## 2006-06-30 0.016739029 -0.142847324 0.085608654
## 2006-07-31 -0.162695331 -0.032124363 -0.231870347
## 2006-08-31 0.209693929 0.145516386 -0.135356116
## 2006-09-29 -0.066146168 0.019601080 -0.018123853
## 2006-10-31 0.195018557 -0.008743114 -0.089378946
## 2006-11-30 0.117824938 0.017150792 0.005238670
## 2006-12-29 -0.083454009 -0.005185211 0.023770224
## 2007-01-31 -0.027300021 -0.048377123 0.201258918
## 2007-02-28 -0.007245408 -0.012121523 -0.009562584
## 2007-03-30 -0.023190702 -0.181669357 0.012521020
## 2007-04-30 0.038291077 -0.019832555 -0.101765593
## 2007-05-31 0.115412615 -0.092835046 0.002103921
## 2007-06-29 0.047441666 0.035574170 -0.198913601
## 2007-07-31 -0.103091238 0.112294669 -0.105876584
## 2007-08-31 -0.013655855 0.119459567 -0.130197299
## 2007-09-28 0.002908976 0.048493022 0.015440312
## 2007-10-31 0.241574545 0.018201445 0.271554940
## 2007-11-30 -0.277306068 0.115020129 0.042469525
## 2007-12-31 -0.095115480 -0.029360378 -0.137711997
## 2008-01-31 -0.405461779 -0.121175344 0.111358251
## 2008-02-29 -0.104296207 -0.073896423 0.145698809
## 2008-03-31 -0.223480518 -0.030966938 0.053367554
## 2008-04-30 0.296666507 0.126638332 0.202377991
## 2008-05-30 0.154807636 -0.031550999 -0.079029445
## 2008-06-30 -0.335756974 0.035039375 0.094487618
## 2008-07-31 -0.285391724 -0.120174799 0.401350806
## 2008-08-29 0.031897972 0.073356926 -0.061116199
## 2008-09-30 -0.413674086 0.059058220 -0.063131782
## 2008-10-31 -0.518795332 0.284102149 0.021651302
## 2008-11-28 -0.532124237 0.074210989 -0.182860353
## 2008-12-31 0.127894855 0.205821985 0.049865952
## 2009-01-30 0.002727966 -0.282904294 -0.132090045
## 2009-02-27 -0.058901878 0.034983327 -0.028166041
## 2009-03-31 0.225219121 0.423092270 0.026185577
## 2009-04-30 0.333561531 0.073838965 0.095017327
## 2009-05-29 0.068055380 -0.079081889 0.147240930
## 2009-06-30 -0.093206280 0.222015263 -0.003657794
## 2009-07-31 0.329838830 -0.074286745 0.120311700
## 2009-08-31 0.051645786 -0.092955621 -0.071085205
## 2009-09-30 0.197442256 -0.064011322 0.065447791
## 2009-10-30 -0.218680655 -0.184804193 -0.023231054
## 2009-11-30 0.099656282 -0.003893188 0.055788533
## 2009-12-31 0.311595167 0.054479549 -0.223711070
# time series plot of principal component factors
# library(PerformanceAnalytics)
# chart.TimeSeries(pc.fit$scores[, 1, drop=FALSE], colorset="blue")
# compare with direct eigen-value analysis
# notice the sign change on the first set of loadings
eigen.fit = eigen(var(X_i))
names(eigen.fit)
## [1] "values" "vectors"
names(eigen.fit$values) = rownames(eigen.fit$vectors) = colnames(X_i)
cbind(pc.fit$loadings[,1:3], eigen.fit$vectors[, 1:3])
## Comp.1 Comp.2 Comp.3
## CERN 0.2432381 0.14806546 0.01738893 -0.2432381 0.14806546 0.01738893
## EQIX 0.4658160 0.16706284 0.20117675 -0.4658160 0.16706284 0.20117675
## HFC 0.3328147 -0.60344658 -0.21841162 -0.3328147 -0.60344658 -0.21841162
## JKHY 0.1232816 0.07003749 0.09042850 -0.1232816 0.07003749 0.09042850
## KSU 0.2882763 -0.50464492 0.33505073 -0.2882763 -0.50464492 0.33505073
## NTES 0.3266215 0.42116958 0.56028310 -0.3266215 0.42116958 0.56028310
## NVDA 0.5367067 0.30963973 -0.67708560 -0.5367067 0.30963973 -0.67708560
## SBAC 0.3464791 -0.23062206 0.13683397 -0.3464791 -0.23062206 0.13683397
compute uncentered pc factors from eigenvectors and return data
pc.factors.uc = X_i %*% eigen.fit$vectors
colnames(pc.factors.uc) = paste(colnames(pc.fit$scores),".uc",sep="")
compare centered and uncentered scores. Note sign change on first factor
We can treat centered scores as unobservable factor values (F)
cbind(pc.fit$scores[,1,drop=F], -pc.factors.uc[, 1, drop=F])
## Comp.1 Comp.1.uc
## 2005-01-31 -0.163950105 -0.102699069
## 2005-02-28 0.233003218 0.294254255
## 2005-03-31 -0.129569915 -0.068318879
## 2005-04-29 -0.212284705 -0.151033669
## 2005-05-31 0.323732402 0.384983438
## 2005-06-30 0.190758905 0.252009942
## 2005-07-29 0.113915948 0.175166985
## 2005-08-31 0.047590584 0.108841620
## 2005-09-30 0.236588793 0.297839829
## 2005-10-31 -0.253968414 -0.192717377
## 2005-11-30 0.119930888 0.181181924
## 2005-12-30 -0.090741775 -0.029490738
## 2006-01-31 0.405708538 0.466959575
## 2006-02-28 0.000738157 0.061989194
## 2006-03-31 0.351376242 0.412627278
## 2006-04-28 -0.083983704 -0.022732667
## 2006-05-31 -0.255789725 -0.194538688
## 2006-06-30 0.016739029 0.077990066
## 2006-07-31 -0.162695331 -0.101444294
## 2006-08-31 0.209693929 0.270944966
## 2006-09-29 -0.066146168 -0.004895131
## 2006-10-31 0.195018557 0.256269593
## 2006-11-30 0.117824938 0.179075975
## 2006-12-29 -0.083454009 -0.022202973
## 2007-01-31 -0.027300021 0.033951015
## 2007-02-28 -0.007245408 0.054005629
## 2007-03-30 -0.023190702 0.038060334
## 2007-04-30 0.038291077 0.099542114
## 2007-05-31 0.115412615 0.176663652
## 2007-06-29 0.047441666 0.108692703
## 2007-07-31 -0.103091238 -0.041840202
## 2007-08-31 -0.013655855 0.047595181
## 2007-09-28 0.002908976 0.064160013
## 2007-10-31 0.241574545 0.302825581
## 2007-11-30 -0.277306068 -0.216055032
## 2007-12-31 -0.095115480 -0.033864443
## 2008-01-31 -0.405461779 -0.344210743
## 2008-02-29 -0.104296207 -0.043045170
## 2008-03-31 -0.223480518 -0.162229481
## 2008-04-30 0.296666507 0.357917543
## 2008-05-30 0.154807636 0.216058673
## 2008-06-30 -0.335756974 -0.274505937
## 2008-07-31 -0.285391724 -0.224140688
## 2008-08-29 0.031897972 0.093149008
## 2008-09-30 -0.413674086 -0.352423050
## 2008-10-31 -0.518795332 -0.457544296
## 2008-11-28 -0.532124237 -0.470873201
## 2008-12-31 0.127894855 0.189145891
## 2009-01-30 0.002727966 0.063979003
## 2009-02-27 -0.058901878 0.002349158
## 2009-03-31 0.225219121 0.286470158
## 2009-04-30 0.333561531 0.394812567
## 2009-05-29 0.068055380 0.129306416
## 2009-06-30 -0.093206280 -0.031955244
## 2009-07-31 0.329838830 0.391089867
## 2009-08-31 0.051645786 0.112896822
## 2009-09-30 0.197442256 0.258693293
## 2009-10-30 -0.218680655 -0.157429619
## 2009-11-30 0.099656282 0.160907318
## 2009-12-31 0.311595167 0.372846203
use first 3 eigen-vectors to compute three factor (with normalization to have pos correlation with market)
note: cannot treat pc as a portfolio b/c weights do not sum to unity
p3 = pc.fit$loadings[, 1:3]
p3 # refers to B
## Comp.1 Comp.2 Comp.3
## CERN 0.2432381 0.14806546 0.01738893
## EQIX 0.4658160 0.16706284 0.20117675
## HFC 0.3328147 -0.60344658 -0.21841162
## JKHY 0.1232816 0.07003749 0.09042850
## KSU 0.2882763 -0.50464492 0.33505073
## NTES 0.3266215 0.42116958 0.56028310
## NVDA 0.5367067 0.30963973 -0.67708560
## SBAC 0.3464791 -0.23062206 0.13683397
colSums(p3)
## Comp.1 Comp.2 Comp.3
## 2.6632340 -0.2227385 0.4456648
create factor mimicking portfolio by normalizing weights to unity
p3 = p3/colSums(p3)
p3
## Comp.1 Comp.2 Comp.3
## CERN 0.09133186 0.33223507 -0.07806883
## EQIX -2.09131374 0.06272931 0.45140825
## HFC 0.74678271 2.70921598 -0.08200992
## JKHY 0.04629018 0.15715286 -0.40598514
## KSU -1.29423682 -0.18948576 0.75179992
## NTES 0.73288614 -1.89087051 0.21037697
## NVDA 0.20152442 0.69478172 3.03982356
## SBAC -1.55554222 -0.08659474 0.30703340
barplot(p3[,1], horiz=F, main="Factor mimicking weights", col="aquamarine3", cex.names = 0.75, las=2)
create first 3 factors
f3 = X_i %*% p3
head(f3)
## Comp.1 Comp.2 Comp.3
## [1,] 0.08512822 0.5922822 -0.18906622
## [2,] -0.03331918 0.8233623 0.91814466
## [3,] 0.06692398 -0.4219316 -0.49491844
## [4,] 0.44859153 -0.2911100 -0.32647122
## [5,] -0.58586909 0.4041709 0.87676800
## [6,] -0.36567626 0.4038035 0.07493234
head(pc.fit$scores[, 1:3])
## Comp.1 Comp.2 Comp.3
## 2005-01-31 -0.1639501 -0.128585797 -0.13853565
## 2005-02-28 0.2330032 -0.117016589 -0.18711190
## 2005-03-31 -0.1295699 -0.005143263 0.17841038
## 2005-04-29 -0.2122847 0.045327324 0.02292385
## 2005-05-31 0.3237324 -0.044139668 -0.08042189
## 2005-06-30 0.1907589 -0.113152694 0.07468680
# chart.TimeSeries(f1, main="First principal component factor", colorset="blue")
# estimate factor betas by multivariate regression
n.obs <- dim(X_i)[1]
X.mat = cbind(rep(1, n.obs), f3)
colnames(X.mat) = c("intercept", "Factor 1", "Factor 2", "Factor 3")
XX.mat = crossprod(X.mat)
multivariate least squares
G.hat = solve(XX.mat)%*%crossprod(X.mat, X_i)
t(G.hat)
## intercept Factor 1 Factor 2 Factor 3
## CERN 0.0132667447 -0.041437681 -0.001681903 0.062413782
## EQIX 0.0002787206 -0.181123478 -0.024559428 0.086667590
## HFC 0.0111770652 -0.006253207 0.250049423 0.012950841
## JKHY 0.0001910669 -0.054936880 -0.009379329 0.015458254
## KSU 0.0026283536 -0.213190080 0.035118684 -0.008117278
## NTES 0.0167457528 0.036348705 -0.174892936 0.148202735
## NVDA -0.0023117073 0.088779605 0.008259739 0.310786935
## SBAC 0.0111791376 -0.200420285 0.042134703 0.010665633
# can also use solve(qr(X.mat), returns.mat)
beta.hat = G.hat[2:4,]
beta.hat
## CERN EQIX HFC JKHY KSU
## Factor 1 -0.041437681 -0.18112348 -0.006253207 -0.054936880 -0.213190080
## Factor 2 -0.001681903 -0.02455943 0.250049423 -0.009379329 0.035118684
## Factor 3 0.062413782 0.08666759 0.012950841 0.015458254 -0.008117278
## NTES NVDA SBAC
## Factor 1 0.0363487 0.088779605 -0.20042028
## Factor 2 -0.1748929 0.008259739 0.04213470
## Factor 3 0.1482027 0.310786935 0.01066563
B
## [,1] [,2] [,3]
## [1,] -0.05361747 0.018608974 -0.002152947
## [2,] -0.10268075 0.020996577 -0.024907970
## [3,] -0.07336302 -0.075841600 0.027041843
## [4,] -0.02717521 0.008802363 -0.011196077
## [5,] -0.06354532 -0.063424136 -0.041483092
## [6,] -0.07199784 0.052932895 -0.069369422
## [7,] -0.11830733 0.038915744 0.083830900
## [8,] -0.07637508 -0.028984746 -0.016941602
E.hat = X_i - X.mat%*%G.hat
diagD.hat = diag(crossprod(E.hat)/(n.obs-4))
compute covariance/correlation matrices with three pc factor
cov.pc3 = t(beta.hat)%*%var(f3)%*%(beta.hat) + diag(diagD.hat)
cov.pc3
## CERN EQIX HFC JKHY KSU
## CERN 0.009170087 0.004347251 0.0018042356 0.0010630480 0.0029948366
## EQIX 0.004347251 0.015323466 0.0022363178 0.0027852962 0.0082430898
## HFC 0.001804236 0.002236318 0.0178692635 0.0003072202 0.0026750652
## JKHY 0.001063048 0.002785296 0.0003072202 0.0045213792 0.0022103079
## KSU 0.002994837 0.008243090 0.0026750652 0.0022103079 0.0154475841
## NTES 0.001842479 0.004041528 -0.0043249290 0.0010098351 0.0009944295
## NVDA 0.005824709 0.010117360 0.0079642528 0.0021509025 0.0048483234
## SBAC 0.003280556 0.008600031 0.0034090876 0.0022641784 0.0075597770
## NTES NVDA SBAC
## CERN 0.0018424788 0.005824709 0.003280556
## EQIX 0.0040415277 0.010117360 0.008600031
## HFC -0.0043249290 0.007964253 0.003409088
## JKHY 0.0010098351 0.002150903 0.002264178
## KSU 0.0009944295 0.004848323 0.007559777
## NTES 0.0176116688 0.006315000 0.001190172
## NVDA 0.0063150002 0.023657445 0.006267871
## SBAC 0.0011901720 0.006267871 0.013075481
Sigma_PCA3
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.008804419 0.005949832 0.0024639851 0.001644973 0.002316192
## [2,] 0.005949832 0.015092330 0.0052669985 0.003254062 0.006226451
## [3,] 0.002463985 0.005266999 0.0175598543 0.001023307 0.008350286
## [4,] 0.001644973 0.003254062 0.0010233073 0.004328908 0.001633023
## [5,] 0.002316192 0.006226451 0.0083502855 0.001633023 0.015041805
## [6,] 0.004994717 0.010232053 -0.0006084137 0.003199156 0.004095561
## [7,] 0.006887038 0.010876926 0.0079948929 0.002618999 0.001572115
## [8,] 0.003592136 0.007655651 0.0073432235 0.002010054 0.007394401
## [,6] [,7] [,8]
## [1,] 0.0049947171 0.006887038 0.003592136
## [2,] 0.0102320532 0.010876926 0.007655651
## [3,] -0.0006084137 0.007994893 0.007343223
## [4,] 0.0031991560 0.002618999 0.002010054
## [5,] 0.0040955609 0.001572115 0.007394401
## [6,] 0.0170421169 0.004762494 0.005139823
## [7,] 0.0047624940 0.023632553 0.006487539
## [8,] 0.0051398229 0.006487539 0.012804919
diag(cov.pc3)
## CERN EQIX HFC JKHY KSU NTES
## 0.009170087 0.015323466 0.017869263 0.004521379 0.015447584 0.017611669
## NVDA SBAC
## 0.023657445 0.013075481
error difference between pca and empirical covariance matrix
norm(cov.pc3 - cov(X_i), "F")
## [1] 0.01525941
norm(Sigma_PCA3 - cov(X_i), "F")
## [1] 0.01146673
==========================================================
Plot Strategy Statistics Side by Side
layout(1:1)
plotbt.strategy.sidebyside(models)
Plot transition maps
layout(1:len(models))
for(m in names(models)) {
plotbt.transition.map(models[[m]]$weight, name=m)
legend('topright', legend = m, bty = 'n')
}
