Final Project
Julie Anne Miranda
6/18/2021
# https://palomar.home.ece.ust.hk/MAFS6010R_lectures/Rsession_factor_models.html#final_comparison_of_covariance_matrix_estimations_via_different_factor_models
# https://systematicinvestor.wordpress.com/2012/03/19/backtesting-asset-allocation-portfolios/
rm(list = ls())
library(lubridate)##
## Attaching package: 'lubridate'## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(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 zoolibrary(lpSolve)#devtools::install_github('joshuaulrich/xts', force = T)
#devtools::install_github('joshuaulrich/quantmod', force = T)
#devtools::install_github('systematicinvestor/SIT.date', force = T)curl::curl_download('https://github.com/systematicinvestor/SIT/raw/master/SIT.tar.gz', 'SIT.tar.gz',mode = 'wb',quiet=T)
install.packages('SIT.tar.gz', repos = NULL, type='source')## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/4.0'
## (as 'lib' is unspecified)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':
##
## closelibrary(PerformanceAnalytics)##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendlibrary(xts)
library(tidyr)5 Credit Services in USA
#Import data
tickers = spl('AXP,CIT,EFX,V,WU')
datas <- new.env()
getSymbols(tickers, src = 'yahoo', from = '2005-01-01', env = datas, auto.assign = T)## [1] "AXP" "CIT" "EFX" "V" "WU"for(i in ls(datas)) datas[[i]] = adjustOHLC(datas[[i]], use.Adjusted=T)
names(datas)## [1] "AXP" ".getSymbols" "CIT" "EFX" "WU"
## [6] "V"datas$EEM## NULLhead(datas$EEM)## NULL#
bt.prep(datas, align='remove.na')
names(datas)## [1] "prices" "AXP" "dates" "weight"
## [5] ".getSymbols" "symbolnames" "CIT" "execution.price"
## [9] "EFX" "WU" "V"#
head(datas$SPY)## NULLhead(datas$prices)## AXP CIT EFX V WU
## 2009-12-10 33.82547 24.26659 26.03877 18.96895 13.31280
## 2009-12-11 34.24587 24.81068 26.17718 18.74312 13.36184
## 2009-12-14 34.70831 23.85643 26.46266 19.53349 13.45993
## 2009-12-15 34.43926 21.93117 26.29829 19.84918 13.29178
## 2009-12-16 34.69991 21.93954 26.34154 20.05426 13.47395
## 2009-12-17 34.01045 21.79724 26.09066 20.05426 13.37555#
prices_monthly <- to.monthly(datas$prices, indexAt = "last", OHLC = FALSE) # indexAt = 'endof'
head(prices_monthly)## AXP CIT EFX V WU
## 2009-12-31 34.06930 23.11144 26.72217 20.15334 13.24903
## 2010-01-29 31.80258 26.63549 27.68241 18.90212 13.03114
## 2010-02-26 32.25015 30.49437 27.94172 19.68055 11.09123
## 2010-03-31 34.99451 32.61215 31.00786 21.00752 11.96482
## 2010-04-30 39.11651 33.98494 29.10236 20.82289 12.87488
## 2010-05-28 33.81559 30.79571 26.23564 16.74688 11.25935monthly.ret <- na.omit(Return.calculate(prices_monthly, method = "discrete"))
head(monthly.ret)## AXP CIT EFX V WU
## 2010-01-29 -0.0665325713 0.15248077 0.035933940 -0.062085184 -0.01644573
## 2010-02-26 0.0140732293 0.14487735 0.009367321 0.041182208 -0.14886700
## 2010-03-31 0.0850962376 0.06944837 0.109733593 0.067425199 0.07876337
## 2010-04-30 0.1177898157 0.04209425 -0.061452157 -0.008788474 0.07606142
## 2010-05-28 -0.1355160333 -0.09384228 -0.098504857 -0.195746739 -0.12547924
## 2010-06-30 0.0002269663 -0.07964125 -0.072396943 -0.023598846 -0.06235511## 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:
## NULLhead(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)## AXP CIT EFX V WU
## 2010-01-29 -0.0665325713 0.15248077 0.035933940 -0.062085184 -0.01644573
## 2010-02-26 0.0140732293 0.14487735 0.009367321 0.041182208 -0.14886700
## 2010-03-31 0.0850962376 0.06944837 0.109733593 0.067425199 0.07876337
## 2010-04-30 0.1177898157 0.04209425 -0.061452157 -0.008788474 0.07606142
## 2010-05-28 -0.1355160333 -0.09384228 -0.098504857 -0.195746739 -0.12547924
## 2010-06-30 0.0002269663 -0.07964125 -0.072396943 -0.023598846 -0.06235511dim(X)## [1] 135 5#
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.67f <- 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.0025dim(f)## [1] 195 3#======================================================================
# Use Professor Palomar's package to compute estimated coefficients
# Based on CAPM model, compute covariance matrix for the 5-asset portfolio by using past 60-month returns from 2010/01 - 2014/12.
#======================================================================
#devtools::install_github("dppalomar/covFactorModel")
library(covFactorModel)
#
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
## AXP 0.016807357 0.12085846
## CIT 0.011965836 -0.08184649
## EFX 0.018907744 0.38720919
## V 0.021244283 0.96969532
## WU 0.004473712 0.13725370#======================================================================
# 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
## AXP 0.016807357 0.12085846
## CIT 0.011965836 -0.08184649
## EFX 0.018907744 0.38720919
## V 0.021244283 0.96969532
## WU 0.004473712 0.13725370# 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## AXP CIT EFX V WU
## AXP 3.272649e-03 -1.080449e-06 5.111517e-06 1.280087e-05 1.811875e-06
## CIT -1.080449e-06 5.417279e-03 -3.461568e-06 -8.668870e-06 -1.227019e-06
## EFX 5.111517e-06 -3.461568e-06 2.711577e-03 4.101173e-05 5.804928e-06
## V 1.280087e-05 -8.668870e-06 4.101173e-05 3.131891e-03 1.453739e-05
## WU 1.811875e-06 -1.227019e-06 5.804928e-06 1.453739e-05 6.325509e-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]
## AXP 0.12085846
## CIT -0.08184649
## EFX 0.38720919
## V 0.96969532
## WU 0.13725370sigeps = 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]
## [1,] 0.003271053 0.000000000 0.0000000 0.000000000 0.000000000
## [2,] 0.000000000 0.005416547 0.0000000 0.000000000 0.000000000
## [3,] 0.000000000 0.000000000 0.0026952 0.000000000 0.000000000
## [4,] 0.000000000 0.000000000 0.0000000 0.003029184 0.000000000
## [5,] 0.000000000 0.000000000 0.0000000 0.000000000 0.006323451cov_1f = sigF*beta_%*%t(beta_)+sigeps
cov_1f # This is covariance matrix computed from single factor model (from lm())## AXP CIT EFX V WU
## AXP 3.272649e-03 -1.080449e-06 5.111517e-06 1.280087e-05 1.811875e-06
## CIT -1.080449e-06 5.417279e-03 -3.461568e-06 -8.668870e-06 -1.227019e-06
## EFX 5.111517e-06 -3.461568e-06 2.711577e-03 4.101173e-05 5.804928e-06
## V 1.280087e-05 -8.668870e-06 4.101173e-05 3.131891e-03 1.453739e-05
## WU 1.811875e-06 -1.227019e-06 5.804928e-06 1.453739e-05 6.325509e-03#=======================================================================
# Backtesting portfolio using SIT package
#=======================================================================
dates <- '2005-01::2021-03'
X <- monthly.ret[dates]
head(X)## AXP CIT EFX V WU
## 2010-01-29 -0.0665325713 0.15248077 0.035933940 -0.062085184 -0.01644573
## 2010-02-26 0.0140732293 0.14487735 0.009367321 0.041182208 -0.14886700
## 2010-03-31 0.0850962376 0.06944837 0.109733593 0.067425199 0.07876337
## 2010-04-30 0.1177898157 0.04209425 -0.061452157 -0.008788474 0.07606142
## 2010-05-28 -0.1355160333 -0.09384228 -0.098504857 -0.195746739 -0.12547924
## 2010-06-30 0.0002269663 -0.07964125 -0.072396943 -0.023598846 -0.06235511dim(X)## [1] 135 5tail(X)## AXP CIT EFX V WU
## 2020-10-30 -0.08607397 0.66290224 -0.12938166 -0.091313594 -0.09286041
## 2020-11-30 0.29975889 0.15025707 0.22459795 0.159357927 0.16049374
## 2020-12-31 0.01956324 0.07196179 0.15542241 0.039838357 -0.01790802
## 2021-01-29 -0.03508679 0.02785518 -0.08157019 -0.116490609 0.01504098
## 2021-02-26 0.16342680 0.23992292 -0.08397974 0.100748883 0.04265827
## 2021-03-31 0.04883288 0.13583241 0.11891523 -0.003107492 0.07202445#
prices_monthly <- to.monthly(datas$prices, indexAt = "last", OHLC = FALSE) # indexAt = 'endof'
head(prices_monthly)## AXP CIT EFX V WU
## 2009-12-31 34.06930 23.11144 26.72217 20.15334 13.24903
## 2010-01-29 31.80258 26.63549 27.68241 18.90212 13.03114
## 2010-02-26 32.25015 30.49437 27.94172 19.68055 11.09123
## 2010-03-31 34.99451 32.61215 31.00786 21.00752 11.96482
## 2010-04-30 39.11651 33.98494 29.10236 20.82289 12.87488
## 2010-05-28 33.81559 30.79571 26.23564 16.74688 11.25935prices_monthly <- prices_monthly[dates]
head(prices_monthly)## AXP CIT EFX V WU
## 2009-12-31 34.06930 23.11144 26.72217 20.15334 13.24903
## 2010-01-29 31.80258 26.63549 27.68241 18.90212 13.03114
## 2010-02-26 32.25015 30.49437 27.94172 19.68055 11.09123
## 2010-03-31 34.99451 32.61215 31.00786 21.00752 11.96482
## 2010-04-30 39.11651 33.98494 29.10236 20.82289 12.87488
## 2010-05-28 33.81559 30.79571 26.23564 16.74688 11.25935tail(prices_monthly)## AXP CIT EFX V WU
## 2020-10-30 90.64872 28.66063 135.7666 180.8926 18.88968
## 2020-11-30 117.82148 32.96709 166.2595 209.7193 21.92135
## 2020-12-31 120.12645 35.33946 192.1000 218.0742 21.52878
## 2021-01-29 115.91160 36.32385 176.4303 192.6706 21.85260
## 2021-02-26 134.85466 45.03877 161.6138 212.0819 22.78479
## 2021-03-31 141.44000 51.15650 180.8321 211.4229 24.42585#*****************************************************************
# Code Strategies
#******************************************************************
# We need to create environment variable required by SIT
# data_m <- new.env()
#
# data_m$EEM <- prices_monthly$EEM
# colnames(data_m$EEM) <- 'Close'
# head(data_m$EEM)
# Change column names to Close which will be consistent with the requirement of SIT package
# data_m$SPY <- prices_monthly$SPY %>% `colnames<-` (c("Close"))
# data_m$QQQ <- prices_monthly$QQQ %>% `colnames<-` (c("Close"))
# data_m$EEM <- prices_monthly$EEM %>% `colnames<-` (c("Close"))
# data_m$IWM <- prices_monthly$IWM %>% `colnames<-` (c("Close"))
# data_m$EFA <- prices_monthly$EFA %>% `colnames<-` (c("Close"))
# data_m$TLT <- prices_monthly$TLT %>% `colnames<-` (c("Close"))
# data_m$IYR <- prices_monthly$IYR %>% `colnames<-` (c("Close"))
# data_m$GLD <- prices_monthly$GLD %>% `colnames<-` (c("Close"))
# names(data_m)
# head(data_m$EEM)
# 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] "WU" "V" "EFX" "CIT" "AXP"# Check if the column name is 'Close'
head(data_m$EEM)## NULL#
bt.prep(data_m, align='remove.na', dates = dates)
names(data_m)## [1] "prices" "execution.price" "weight" "dates"
## [5] "symbolnames" "WU" "V" "EFX"
## [9] "CIT" "AXP"#===============================================================
# 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)## AXP CIT EFX V WU
## 2009-12-31 0.2 0.2 0.2 0.2 0.2
## 2010-01-29 0.2 0.2 0.2 0.2 0.2
## 2010-02-26 0.2 0.2 0.2 0.2 0.2
## 2010-03-31 0.2 0.2 0.2 0.2 0.2
## 2010-04-30 0.2 0.2 0.2 0.2 0.2
## 2010-05-28 0.2 0.2 0.2 0.2 0.2data_m$weight[1:59, ] <- NA
models$equal.weight = bt.run.share(data_m)## Latest weights :
## AXP CIT EFX V WU
## 2021-03-31 16.35 11.01 22.13 34.72 15.79
##
## Performance summary :
## CAGR Best Worst
## 6.8 19.3 -19.9 # 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 :
## AXP CIT EFX V WU
## 2021-03-31 20 20 20 20 20
##
## Performance summary :
## CAGR Best Worst
## 7.6 19.9 -24.8 # 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':
##
## crossdim(weight)## [1] 136 5# 195 8
head(weight, 70)## AXP CIT EFX V WU
## [1,] NA NA NA NA NA
## [2,] NA NA NA NA NA
## [3,] NA NA NA NA NA
## [4,] NA NA NA NA NA
## [5,] NA NA NA NA NA
## [6,] NA NA NA NA NA
## [7,] NA NA NA NA NA
## [8,] NA NA NA NA NA
## [9,] NA NA NA NA NA
## [10,] NA NA NA NA NA
## [11,] NA NA NA NA NA
## [12,] NA NA NA NA NA
## [13,] NA NA NA NA NA
## [14,] NA NA NA NA NA
## [15,] NA NA NA NA NA
## [16,] NA NA NA NA NA
## [17,] NA NA NA NA NA
## [18,] NA NA NA NA NA
## [19,] NA NA NA NA NA
## [20,] NA NA NA NA NA
## [21,] NA NA NA NA NA
## [22,] NA NA NA NA NA
## [23,] NA NA NA NA NA
## [24,] NA NA NA NA NA
## [25,] NA NA NA NA NA
## [26,] NA NA NA NA NA
## [27,] NA NA NA NA NA
## [28,] NA NA NA NA NA
## [29,] NA NA NA NA NA
## [30,] NA NA NA NA NA
## [31,] NA NA NA NA NA
## [32,] NA NA NA NA NA
## [33,] NA NA NA NA NA
## [34,] NA NA NA NA NA
## [35,] NA NA NA NA NA
## [36,] NA NA NA NA NA
## [37,] NA NA NA NA NA
## [38,] NA NA NA NA NA
## [39,] NA NA NA NA NA
## [40,] NA NA NA NA NA
## [41,] NA NA NA NA NA
## [42,] NA NA NA NA NA
## [43,] NA NA NA NA NA
## [44,] NA NA NA NA NA
## [45,] NA NA NA NA NA
## [46,] NA NA NA NA NA
## [47,] NA NA NA NA NA
## [48,] NA NA NA NA NA
## [49,] NA NA NA NA NA
## [50,] NA NA NA NA NA
## [51,] NA NA NA NA NA
## [52,] NA NA NA NA NA
## [53,] NA NA NA NA NA
## [54,] NA NA NA NA NA
## [55,] NA NA NA NA NA
## [56,] NA NA NA NA NA
## [57,] NA NA NA NA NA
## [58,] NA NA NA NA NA
## [59,] NA NA NA NA NA
## [60,] 0.1736106 0.059873006 0.4248806 0.2674951 0.07414069
## [61,] 0.1738296 0.059154950 0.4256130 0.2678023 0.07360016
## [62,] 0.1526054 0.029719914 0.4453575 0.2877056 0.08461152
## [63,] 0.1937821 0.032511268 0.4284116 0.2844726 0.06082250
## [64,] 0.1918677 0.021607029 0.4490965 0.2602499 0.07717874
## [65,] 0.1716463 0.019729115 0.4627220 0.2702979 0.07560465
## [66,] 0.1446470 -0.026955948 0.4149326 0.3895143 0.07786202
## [67,] 0.1165741 0.002779974 0.4287493 0.3842177 0.06767894
## [68,] 0.1910198 -0.001956280 0.4482710 0.2927450 0.06992053
## [69,] 0.2772466 -0.056848356 0.4532497 0.2681181 0.05823397
## [70,] 0.2650579 -0.057946862 0.4500397 0.2670508 0.07579839tail(weight)## AXP CIT EFX V WU
## [131,] -0.09608713 0.03318857 -0.07768133 0.7412319 0.3993480
## [132,] -0.24161637 0.05802012 -0.12005845 0.8681708 0.4354839
## [133,] -0.24457500 0.06087603 -0.12090505 0.8685630 0.4360410
## [134,] -0.21318118 0.06302270 -0.11748719 0.7901915 0.4774542
## [135,] -0.26366508 0.05792796 -0.07493431 0.7853584 0.4953130
## [136,] -0.24980958 0.05491734 -0.07162575 0.7741723 0.4923457# 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 :
## AXP CIT EFX V WU
## 2021-03-31 -26.37 5.79 -7.49 78.54 49.53
##
## Performance summary :
## CAGR Best Worst
## 6 14.1 -25.5 #=================================================================
# 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] 136 5# 195 8
head(weight_1, 70)## AXP CIT EFX V WU
## [1,] NA NA NA NA NA
## [2,] NA NA NA NA NA
## [3,] NA NA NA NA NA
## [4,] NA NA NA NA NA
## [5,] NA NA NA NA NA
## [6,] NA NA NA NA NA
## [7,] NA NA NA NA NA
## [8,] NA NA NA NA NA
## [9,] NA NA NA NA NA
## [10,] NA NA NA NA NA
## [11,] NA NA NA NA NA
## [12,] NA NA NA NA NA
## [13,] NA NA NA NA NA
## [14,] NA NA NA NA NA
## [15,] NA NA NA NA NA
## [16,] NA NA NA NA NA
## [17,] NA NA NA NA NA
## [18,] NA NA NA NA NA
## [19,] NA NA NA NA NA
## [20,] NA NA NA NA NA
## [21,] NA NA NA NA NA
## [22,] NA NA NA NA NA
## [23,] NA NA NA NA NA
## [24,] NA NA NA NA NA
## [25,] NA NA NA NA NA
## [26,] NA NA NA NA NA
## [27,] NA NA NA NA NA
## [28,] NA NA NA NA NA
## [29,] NA NA NA NA NA
## [30,] NA NA NA NA NA
## [31,] NA NA NA NA NA
## [32,] NA NA NA NA NA
## [33,] NA NA NA NA NA
## [34,] NA NA NA NA NA
## [35,] NA NA NA NA NA
## [36,] NA NA NA NA NA
## [37,] NA NA NA NA NA
## [38,] NA NA NA NA NA
## [39,] NA NA NA NA NA
## [40,] NA NA NA NA NA
## [41,] NA NA NA NA NA
## [42,] NA NA NA NA NA
## [43,] NA NA NA NA NA
## [44,] NA NA NA NA NA
## [45,] NA NA NA NA NA
## [46,] NA NA NA NA NA
## [47,] NA NA NA NA NA
## [48,] NA NA NA NA NA
## [49,] NA NA NA NA NA
## [50,] NA NA NA NA NA
## [51,] NA NA NA NA NA
## [52,] NA NA NA NA NA
## [53,] NA NA NA NA NA
## [54,] NA NA NA NA NA
## [55,] NA NA NA NA NA
## [56,] NA NA NA NA NA
## [57,] NA NA NA NA NA
## [58,] NA NA NA NA NA
## [59,] NA NA NA NA NA
## [60,] 0.2267588 0.1397253 0.2780747 0.2378608 0.11758039
## [61,] 0.2268346 0.1392522 0.2785739 0.2381457 0.11719361
## [62,] 0.2122614 0.1447678 0.2792257 0.2463924 0.11735277
## [63,] 0.2134362 0.1536996 0.2681911 0.2452820 0.11939110
## [64,] 0.2124237 0.1526921 0.2769904 0.2394656 0.11842828
## [65,] 0.2184923 0.1493418 0.2818420 0.2341458 0.11617817
## [66,] 0.2154971 0.1339997 0.2692085 0.2765079 0.10478686
## [67,] 0.2118094 0.1356342 0.2771225 0.2725028 0.10293116
## [68,] 0.2233263 0.1350599 0.2890457 0.2500243 0.10254383
## [69,] 0.2387502 0.1287443 0.2891489 0.2454133 0.09794334
## [70,] 0.2376034 0.1299242 0.2866625 0.2444721 0.10133773tail(weight_1)## AXP CIT EFX V WU
## [131,] 0.1965798 0.04556283 0.1435287 0.3495973 0.2647313
## [132,] 0.1744289 0.05226084 0.1469375 0.3554550 0.2709178
## [133,] 0.1747856 0.05199667 0.1397882 0.3574815 0.2759481
## [134,] 0.2057173 0.05450894 0.1367237 0.3262764 0.2767737
## [135,] 0.1978111 0.05377082 0.1378753 0.3276153 0.2829275
## [136,] 0.2015658 0.05322949 0.1362282 0.3281289 0.2808477# 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 :
## AXP CIT EFX V WU
## 2021-03-31 19.78 5.38 13.79 32.76 28.29
##
## Performance summary :
## CAGR Best Worst
## 7.4 19.6 -20.5 models$mvp.capm.cov$cagr## [1] 0.0737516models$mvp.hist.cov$cagr## [1] 0.05963121#=================================================================
# 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] 136 5# 195 8
head(weight_3, 70)## AXP CIT EFX V WU
## [1,] NA NA NA NA NA
## [2,] NA NA NA NA NA
## [3,] NA NA NA NA NA
## [4,] NA NA NA NA NA
## [5,] NA NA NA NA NA
## [6,] NA NA NA NA NA
## [7,] NA NA NA NA NA
## [8,] NA NA NA NA NA
## [9,] NA NA NA NA NA
## [10,] NA NA NA NA NA
## [11,] NA NA NA NA NA
## [12,] NA NA NA NA NA
## [13,] NA NA NA NA NA
## [14,] NA NA NA NA NA
## [15,] NA NA NA NA NA
## [16,] NA NA NA NA NA
## [17,] NA NA NA NA NA
## [18,] NA NA NA NA NA
## [19,] NA NA NA NA NA
## [20,] NA NA NA NA NA
## [21,] NA NA NA NA NA
## [22,] NA NA NA NA NA
## [23,] NA NA NA NA NA
## [24,] NA NA NA NA NA
## [25,] NA NA NA NA NA
## [26,] NA NA NA NA NA
## [27,] NA NA NA NA NA
## [28,] NA NA NA NA NA
## [29,] NA NA NA NA NA
## [30,] NA NA NA NA NA
## [31,] NA NA NA NA NA
## [32,] NA NA NA NA NA
## [33,] NA NA NA NA NA
## [34,] NA NA NA NA NA
## [35,] NA NA NA NA NA
## [36,] NA NA NA NA NA
## [37,] NA NA NA NA NA
## [38,] NA NA NA NA NA
## [39,] NA NA NA NA NA
## [40,] NA NA NA NA NA
## [41,] NA NA NA NA NA
## [42,] NA NA NA NA NA
## [43,] NA NA NA NA NA
## [44,] NA NA NA NA NA
## [45,] NA NA NA NA NA
## [46,] NA NA NA NA NA
## [47,] NA NA NA NA NA
## [48,] NA NA NA NA NA
## [49,] NA NA NA NA NA
## [50,] NA NA NA NA NA
## [51,] NA NA NA NA NA
## [52,] NA NA NA NA NA
## [53,] NA NA NA NA NA
## [54,] NA NA NA NA NA
## [55,] NA NA NA NA NA
## [56,] NA NA NA NA NA
## [57,] NA NA NA NA NA
## [58,] NA NA NA NA NA
## [59,] NA NA NA NA NA
## [60,] 0.2268474 0.1345365 0.2845171 0.2394584 0.11464066
## [61,] 0.2268753 0.1337620 0.2853501 0.2398765 0.11413613
## [62,] 0.2103763 0.1398293 0.2874491 0.2482942 0.11405107
## [63,] 0.2105597 0.1500175 0.2761357 0.2473076 0.11597953
## [64,] 0.2071415 0.1477333 0.2865434 0.2413339 0.11724785
## [65,] 0.2153305 0.1449013 0.2890256 0.2347371 0.11600544
## [66,] 0.2095910 0.1283711 0.2795525 0.2784572 0.10402826
## [67,] 0.2059771 0.1294575 0.2884884 0.2751013 0.10097571
## [68,] 0.2183147 0.1291350 0.3028389 0.2490450 0.10066639
## [69,] 0.2355544 0.1221977 0.3014564 0.2441741 0.09661736
## [70,] 0.2334938 0.1221313 0.3011809 0.2437870 0.09940704tail(weight_3)## AXP CIT EFX V WU
## [131,] 0.1949617 0.04995576 0.1386755 0.3504222 0.2659848
## [132,] 0.1738424 0.05668030 0.1431061 0.3551733 0.2711979
## [133,] 0.1740234 0.05646848 0.1362283 0.3575157 0.2757642
## [134,] 0.2042013 0.06201688 0.1304179 0.3241251 0.2792388
## [135,] 0.1959575 0.06001338 0.1337282 0.3252141 0.2850868
## [136,] 0.1992377 0.05980865 0.1324103 0.3256792 0.2828642# 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 :
## AXP CIT EFX V WU
## 2021-03-31 19.6 6 13.37 32.52 28.51
##
## Performance summary :
## CAGR Best Worst
## 7.4 19.6 -20.6 #============================================================
# 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
## AXP 0.016732626 -0.03724781 -0.000884161 -0.02310065
## CIT 0.012016445 -0.05820076 -0.035763338 0.02288724
## EFX 0.018668319 -0.02408505 -0.011065154 -0.01171002
## V 0.020644689 -0.03045785 -0.007760490 -0.03044779
## WU 0.004388844 -0.06683948 0.039157384 0.01103845#
# 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.04182299 0.013416021 -0.0224280442
## [2,] -0.05720970 0.009805273 0.0429559985
## [3,] -0.02859771 0.024966469 0.0004055309
## [4,] -0.03413242 0.026554585 -0.0230637764
## [5,] -0.06505482 -0.042155380 -0.0114344636# fiti = lm(formula = X_i ~ t(B))
Sigma_PCA3 <- Sigma
Sigma_PCA3## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.003217207 0.001560810 0.0015218970 0.002301052 0.0024116822
## [2,] 0.001560810 0.005325473 0.0018982895 0.001222353 0.0028172433
## [3,] 0.001521897 0.001898289 0.0026658952 0.001629730 0.0008033108
## [4,] 0.002301052 0.001222353 0.0016297300 0.003080549 0.0013647815
## [5,] 0.002411682 0.002817243 0.0008033108 0.001364781 0.0062183316diag(Sigma_PCA3)## [1] 0.003217207 0.005325473 0.002665895 0.003080549 0.006218332# error
norm(Sigma_PCA3 - cov(X_i), "F")## [1] 0.00127828#---------------------------------------------------------------
# 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
## 0.10511807 0.05766557 0.05441415 0.04026810 0.03475800
##
## 5 variables and 60 observations.summary(pc.fit)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 0.1051181 0.05766557 0.05441415 0.04026810 0.0347580
## Proportion of Variance 0.5479516 0.16490001 0.14682876 0.08040993 0.0599097
## Cumulative Proportion 0.5479516 0.71285161 0.85968037 0.94009030 1.0000000#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
## 0.10511807 0.05766557 0.05441415 0.04026810 0.03475800
##
## 5 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.04182299 0.013416021 -0.0224280442
## [2,] -0.05720970 0.009805273 0.0429559985
## [3,] -0.02859771 0.024966469 0.0004055309
## [4,] -0.03413242 0.026554585 -0.0230637764
## [5,] -0.06505482 -0.042155380 -0.0114344636#
#
plot(pc.fit)
loadings(pc.fit)##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## AXP 0.395 0.231 0.409 0.209 0.762
## CIT 0.540 0.169 -0.783 0.259
## EFX 0.270 0.429 -0.861
## V 0.322 0.457 0.420 0.346 -0.625
## WU 0.614 -0.725 0.208 -0.165 -0.165
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## SS loadings 1.0 1.0 1.0 1.0 1.0
## Proportion Var 0.2 0.2 0.2 0.2 0.2
## Cumulative Var 0.2 0.4 0.6 0.8 1.0pc.fit$loadings[, 1:3]## Comp.1 Comp.2 Comp.3
## AXP 0.3945373 0.2307053 0.408723827
## CIT 0.5396879 0.1686139 -0.782820825
## EFX 0.2697766 0.4293297 -0.007390308
## V 0.3219883 0.4566393 0.420309273
## WU 0.6136949 -0.7249145 0.208379191# 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
## 2010-01-29 0.008189417 -0.01078726 -0.18323209
## 2010-02-26 -0.019294482 0.13827109 -0.12832772
## 2010-03-31 0.143240626 0.03199934 0.01747032
## 2010-04-30 0.068996788 -0.07140896 0.02091511
## 2010-05-28 -0.298183908 -0.10794932 -0.09650642
## 2010-06-30 -0.135752250 -0.03017940 0.03317428pc.fit$scores[, 1:3]## Comp.1 Comp.2 Comp.3
## 2010-01-29 0.008189417 -0.010787261 -0.1832320907
## 2010-02-26 -0.019294482 0.138271087 -0.1283277216
## 2010-03-31 0.143240626 0.031999340 0.0174703185
## 2010-04-30 0.068996788 -0.071408961 0.0209151075
## 2010-05-28 -0.298183908 -0.107949317 -0.0965064206
## 2010-06-30 -0.135752250 -0.030179398 0.0331742764
## 2010-07-30 0.159175706 0.023798921 0.0192047991
## 2010-08-31 -0.120108648 -0.070321090 -0.0884626853
## 2010-09-30 0.175586467 -0.023288929 -0.0141023041
## 2010-10-29 0.033029570 0.041989491 -0.0379859699
## 2010-11-30 -0.060966173 -0.030539021 0.0582272354
## 2010-12-31 0.101677971 -0.039387051 -0.1694317227
## 2011-01-31 0.040161313 -0.083348003 0.0053855311
## 2011-02-28 -0.006588653 -0.071562176 0.1043096474
## 2011-03-31 -0.028480156 0.064183845 0.0140622696
## 2011-04-29 0.031664546 -0.003131279 0.0623785269
## 2011-05-31 0.012426989 0.045587112 -0.0105952653
## 2011-06-30 -0.048754991 -0.020112937 0.0092591169
## 2011-07-29 -0.112149073 -0.019849875 0.0599467519
## 2011-08-31 -0.196802082 0.054505710 0.0738861197
## 2011-09-30 -0.195235296 -0.044683333 0.0241567637
## 2011-10-31 0.258784044 0.033601774 -0.0027558128
## 2011-11-30 -0.033711481 0.010330197 0.0114342168
## 2011-12-30 0.040004877 -0.017042988 -0.0073828337
## 2012-01-31 0.076564609 -0.023694804 -0.0470820664
## 2012-02-29 0.051693168 0.173901213 0.0110064971
## 2012-03-30 0.043392551 0.023345353 0.0295430764
## 2012-04-30 -0.004087038 -0.021529467 0.1022483824
## 2012-05-31 -0.197020705 -0.007612335 -0.0091397205
## 2012-06-29 0.064708426 0.019793747 0.0147748819
## 2012-07-31 0.020860929 -0.020281501 -0.0025915886
## 2012-08-31 -0.005352992 -0.029768829 -0.0289580372
## 2012-09-28 0.030614999 -0.018739216 -0.0227627770
## 2012-10-31 -0.216142737 0.233305344 -0.0180602286
## 2012-11-30 -0.002528713 0.034682708 0.0282307034
## 2012-12-31 0.081018351 -0.041505320 -0.0054234292
## 2013-01-31 0.099125874 0.024583448 -0.0447119115
## 2013-02-28 -0.033205942 -0.020316394 0.0257124961
## 2013-03-28 0.111852024 -0.001376150 0.0440252331
## 2013-04-30 -0.027851644 0.014322108 0.0104981665
## 2013-05-31 0.143767007 -0.031690936 0.0180538317
## 2013-06-28 0.005829391 -0.061131635 0.0006901823
## 2013-07-31 0.048835128 -0.028951930 -0.0728792837
## 2013-08-30 -0.097130224 -0.048692692 0.0078866948
## 2013-09-30 0.080659652 -0.008070169 0.0539234661
## 2013-10-31 -0.022932709 0.109285384 0.0315799839
## 2013-11-29 0.030184296 0.050381527 -0.0155235666
## 2013-12-31 0.076439312 0.022121293 0.0393936525
## 2014-01-31 -0.181620298 0.016747903 0.0160367413
## 2014-02-28 0.097756503 -0.032994821 0.0249360306
## 2014-03-31 -0.060311922 -0.044270685 -0.0395551004
## 2014-04-30 -0.130736975 -0.035551204 0.0453088010
## 2014-05-30 0.042726338 0.012791715 0.0143114151
## 2014-06-30 0.052566525 -0.063382861 -0.0058034734
## 2014-07-31 0.002745927 -0.007693564 -0.0909484661
## 2014-08-29 -0.018238218 0.001111672 0.0198102936
## 2014-09-30 -0.117420480 0.002027059 0.0029853398
## 2014-10-31 0.100633887 0.022024613 0.0221572957
## 2014-11-28 0.079703203 -0.027304925 0.0535697699
## 2014-12-31 -0.044008626 0.013458497 0.0117288596# 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
## AXP 0.3945373 0.2307053 0.408723827 -0.3945373 0.2307053 -0.408723827
## CIT 0.5396879 0.1686139 -0.782820825 -0.5396879 0.1686139 0.782820825
## EFX 0.2697766 0.4293297 -0.007390308 -0.2697766 0.4293297 0.007390308
## V 0.3219883 0.4566393 0.420309273 -0.3219883 0.4566393 -0.420309273
## WU 0.6136949 -0.7249145 0.208379191 -0.6136949 -0.7249145 -0.208379191# 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
## 2010-01-29 0.008189417 0.0356532261
## 2010-02-26 -0.019294482 0.0081693275
## 2010-03-31 0.143240626 0.1707044352
## 2010-04-30 0.068996788 0.0964605973
## 2010-05-28 -0.298183908 -0.2707200988
## 2010-06-30 -0.135752250 -0.1082884413
## 2010-07-30 0.159175706 0.1866395149
## 2010-08-31 -0.120108648 -0.0926448391
## 2010-09-30 0.175586467 0.2030502759
## 2010-10-29 0.033029570 0.0604933794
## 2010-11-30 -0.060966173 -0.0335023641
## 2010-12-31 0.101677971 0.1291417798
## 2011-01-31 0.040161313 0.0676251222
## 2011-02-28 -0.006588653 0.0208751561
## 2011-03-31 -0.028480156 -0.0010163467
## 2011-04-29 0.031664546 0.0591283552
## 2011-05-31 0.012426989 0.0398907980
## 2011-06-30 -0.048754991 -0.0212911818
## 2011-07-29 -0.112149073 -0.0846852642
## 2011-08-31 -0.196802082 -0.1693382735
## 2011-09-30 -0.195235296 -0.1677714867
## 2011-10-31 0.258784044 0.2862478533
## 2011-11-30 -0.033711481 -0.0062476719
## 2011-12-30 0.040004877 0.0674686859
## 2012-01-31 0.076564609 0.1040284179
## 2012-02-29 0.051693168 0.0791569770
## 2012-03-30 0.043392551 0.0708563600
## 2012-04-30 -0.004087038 0.0233767713
## 2012-05-31 -0.197020705 -0.1695568957
## 2012-06-29 0.064708426 0.0921722350
## 2012-07-31 0.020860929 0.0483247377
## 2012-08-31 -0.005352992 0.0221108169
## 2012-09-28 0.030614999 0.0580788076
## 2012-10-31 -0.216142737 -0.1886789283
## 2012-11-30 -0.002528713 0.0249350959
## 2012-12-31 0.081018351 0.1084821602
## 2013-01-31 0.099125874 0.1265896828
## 2013-02-28 -0.033205942 -0.0057421333
## 2013-03-28 0.111852024 0.1393158333
## 2013-04-30 -0.027851644 -0.0003878347
## 2013-05-31 0.143767007 0.1712308165
## 2013-06-28 0.005829391 0.0332931999
## 2013-07-31 0.048835128 0.0762989375
## 2013-08-30 -0.097130224 -0.0696664147
## 2013-09-30 0.080659652 0.1081234608
## 2013-10-31 -0.022932709 0.0045311002
## 2013-11-29 0.030184296 0.0576481055
## 2013-12-31 0.076439312 0.1039031214
## 2014-01-31 -0.181620298 -0.1541564893
## 2014-02-28 0.097756503 0.1252203125
## 2014-03-31 -0.060311922 -0.0328481129
## 2014-04-30 -0.130736975 -0.1032731658
## 2014-05-30 0.042726338 0.0701901473
## 2014-06-30 0.052566525 0.0800303341
## 2014-07-31 0.002745927 0.0302097359
## 2014-08-29 -0.018238218 0.0092255907
## 2014-09-30 -0.117420480 -0.0899566708
## 2014-10-31 0.100633887 0.1280976960
## 2014-11-28 0.079703203 0.1071670118
## 2014-12-31 -0.044008626 -0.0165448169#
# use first 3 eigen-vectors to compue 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
## AXP 0.3945373 0.2307053 0.408723827
## CIT 0.5396879 0.1686139 -0.782820825
## EFX 0.2697766 0.4293297 -0.007390308
## V 0.3219883 0.4566393 0.420309273
## WU 0.6136949 -0.7249145 0.208379191colSums(p3)## Comp.1 Comp.2 Comp.3
## 2.1396850 0.5603736 0.2472012# create factor mimicking portfolio by normalizing weights to unity
p3 = p3/colSums(p3)
p3## Comp.1 Comp.2 Comp.3
## AXP 0.1843904 0.93326932 0.729377314
## CIT 0.9630859 0.07880316 -3.166736086
## EFX 1.0913241 0.76614891 -0.003453924
## V 0.1504840 1.84723771 0.750051816
## WU 1.0951530 -0.33879497 0.842953948barplot(p3[,1], horiz=F, main="Factor mimicking weights", col="blue", cex.names = 0.75, las=2)
# create first 3 factors
f3 = X_i %*% p3
head(f3)## Comp.1 Comp.2 Comp.3
## [1,] 0.146446268 -0.13166045 -0.59194791
## [2,] -0.004487801 0.15823639 -0.54315522
## [3,] 0.298734927 0.26682847 -0.04127011
## [4,] 0.077171825 0.02416165 0.01034856
## [5,] -0.389742459 -0.52841644 -0.05392176
## [6,] -0.227507680 -0.08399808 0.18235557head(pc.fit$scores[, 1:3])## Comp.1 Comp.2 Comp.3
## 2010-01-29 0.008189417 -0.01078726 -0.18323209
## 2010-02-26 -0.019294482 0.13827109 -0.12832772
## 2010-03-31 0.143240626 0.03199934 0.01747032
## 2010-04-30 0.068996788 -0.07140896 0.02091511
## 2010-05-28 -0.298183908 -0.10794932 -0.09650642
## 2010-06-30 -0.135752250 -0.03017940 0.03317428# 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
## AXP -0.001500268 0.13002109 0.19368183 0.063476179
## CIT -0.001523326 0.18798351 0.05685706 -0.241979512
## EFX 0.004749503 0.12217519 0.12782575 -0.006238398
## V -0.001680584 -0.01324126 0.34094848 0.018580885
## WU -0.002909746 0.60598014 -0.25683890 0.205774574# can also use solve(qr(X.mat), returns.mat)
beta.hat = G.hat[2:4,]
beta.hat## AXP CIT EFX V WU
## Factor 1 0.13002109 0.18798351 0.122175193 -0.01324126 0.6059801
## Factor 2 0.19368183 0.05685706 0.127825753 0.34094848 -0.2568389
## Factor 3 0.06347618 -0.24197951 -0.006238398 0.01858088 0.2057746B## [,1] [,2] [,3]
## [1,] -0.04182299 0.013416021 -0.0224280442
## [2,] -0.05720970 0.009805273 0.0429559985
## [3,] -0.02859771 0.024966469 0.0004055309
## [4,] -0.03413242 0.026554585 -0.0230637764
## [5,] -0.06505482 -0.042155380 -0.0114344636E.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## AXP CIT EFX V WU
## AXP 0.003275959 0.001473543 0.001610502 0.002206241 0.002093017
## CIT 0.001473543 0.005332755 0.001940727 0.001244919 0.002500599
## EFX 0.001610502 0.001940727 0.002734488 0.001609023 0.001827145
## V 0.002206241 0.001244919 0.001609023 0.003103769 0.001045026
## WU 0.002093017 0.002500599 0.001827145 0.001045026 0.006245163Sigma_PCA3## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.003217207 0.001560810 0.0015218970 0.002301052 0.0024116822
## [2,] 0.001560810 0.005325473 0.0018982895 0.001222353 0.0028172433
## [3,] 0.001521897 0.001898289 0.0026658952 0.001629730 0.0008033108
## [4,] 0.002301052 0.001222353 0.0016297300 0.003080549 0.0013647815
## [5,] 0.002411682 0.002817243 0.0008033108 0.001364781 0.0062183316diag(cov.pc3)## AXP CIT EFX V WU
## 0.003275959 0.005332755 0.002734488 0.003103769 0.006245163# error difference between pca and empirical covariance matrix
norm(cov.pc3 - cov(X_i), "F")## [1] 0.001755478norm(Sigma_PCA3 - cov(X_i), "F")## [1] 0.00127828#==========================================================
# Plot perfromance
plotbt(models, plotX = T, log = 'y', LeftMargin = 3)
mtext('Cumulative Performance', side = 2, line = 1)
# 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')
}
#