CoCalc -- demo_robustCovMatForPA.R

GitHub Repository: braverock/portfolioanalytics
Path: blob/master/demo/demo_robustCovMatForPA.R
¹⁴³³ views
1
## ----setup, include=FALSE-------------------------------
2
knitr::opts_chunk$set(echo = TRUE, fig.width = 6, fig.height = 4, fig.align = "center")
3

4

5
## ---- warning = FALSE, message = FALSE------------------
6
library(PCRA)
7
library(PortfolioAnalytics)
8
library(CVXR)
9
library(xts)
10

11

12
## -------------------------------------------------------
13
stocksCRSPweekly <- getPCRAData("stocksCRSPweekly")
14
dateRange    <- c("2006-01-01", "2012-12-31")
15
stockItems <- c("Date", "TickerLast", "CapGroupLast", "Return", "MktIndexCRSP")
16
returnsAll <- selectCRSPandSPGMI("weekly",
17
                                   dateRange = dateRange,
18
                                   stockItems = stockItems, 
19
                                   factorItems = NULL, 
20
                                   subsetType = "CapGroupLast",
21
                                   subsetValues = "SmallCap", 
22
                                   outputType = "xts")
23
returns <- returnsAll[,1:30]
24
MARKET <- returnsAll[, 107]
25
returns10 <- returnsAll[140:300, 21:30]
26
range(index(returns))
27
range(index(returns10))
28

29

30
## -------------------------------------------------------
31
tsPlotMP(MARKET)
32

33

34
## -------------------------------------------------------
35
tsPlotMP(returns[, 21:30])
36

37

38
## -------------------------------------------------------
39
tsPlotMP(returns10)
40

41

42
## -------------------------------------------------------
43
funds <- colnames(returns10)
44
pspec <- portfolio.spec(funds)
45
pspec <- add.constraint(pspec, type="full_investment")
46
pspec <- add.constraint(pspec, type="long_only")
47
pspec <- add.objective(pspec, type="risk", name="var")
48

49

50
## -------------------------------------------------------
51
opt <- optimize.portfolio(returns10, pspec, optimize_method = "CVXR")
52

53

54
## -------------------------------------------------------
55
opt
56

57

58
## -------------------------------------------------------
59
outCovClassic <- opt.outputMvo(opt, returns10, digits = 3, frequency = "weekly")
60
(WtsCovClassic <- outCovClassic$Wgts)
61

62

63
## ---- eval = FALSE--------------------------------------
64
## custom.covRob.MM <- function(R, ...){
65
##   out <- list()
66
##   if(hasArg(tol)) tol = match.call(expand.dots = TRUE)$tol else tol = 1e-4
67
##   if(hasArg(maxit)) maxit = match.call(expand.dots = TRUE)$maxit else maxit = 50
68
##   robustCov <- RobStatTM::covRobMM(X=R, tol=tol, maxit=maxit)
69
##   out$sigma <- robustCov$cov
70
##   out$mu <- robustCov$center
71
##   return(out)
72
## }
73

74

75
## -------------------------------------------------------
76
opt <- optimize.portfolio(returns10, pspec, 
77
                   optimize_method = "CVXR", 
78
                   momentFUN = "custom.covRob.MM",
79
                   maxit = 100, tol = 1e-5)
80
outCovRobMM <- opt.outputMvo(opt, returns10, digits = 3, frequency = "weekly")
81
(WtsCovRobMM <- outCovRobMM$Wgts)
82

83

84
## ---- eval = FALSE--------------------------------------
85
## custom.covRob.Rocke <- function(R, ...){
86
##   out <- list()
87
##   if(hasArg(tol)) tol = match.call(expand.dots = TRUE)$tol else tol = 1e-4
88
##   if(hasArg(maxit)) maxit = match.call(expand.dots = TRUE)$maxit else maxit = 50
89
##   if(hasArg(initial)) initial = match.call(expand.dots = TRUE)$initial else initial = 'K'
90
##   if(hasArg(maxsteps)) maxsteps = match.call(expand.dots = TRUE)$maxsteps else maxsteps = 5
91
##   if(hasArg(propmin)) propmin = match.call(expand.dots = TRUE)$propmin else propmin = 2
92
##   if(hasArg(qs)) qs = match.call(expand.dots = TRUE)$qs else qs = 50
93
## 
94
##   robustCov <- RobStatTM::covRobRocke(X = R, initial = initial, maxsteps = maxsteps, propmin = propmin,
95
##                                       qs = qs, tol = tol, maxit = maxit)
96
## 
97
##   out$sigma <- robustCov$cov
98
##   out$mu <- robustCov$center
99
##   return(out)
100
## }
101

102

103
## -------------------------------------------------------
104
opt <- optimize.portfolio(returns10, pspec, 
105
                   optimize_method = "CVXR", 
106
                   momentFUN = "custom.covRob.Rocke",
107
                   tol = 1e-5, maxit =100, maxsteps = 7)
108
outCovRobRocke <- opt.outputMvo(opt, returns10, digits = 3, frequency = "weekly")
109
(WtsCovRobRocke <- outCovRobRocke$Wgts)
110

111

112
## ---- eval = FALSE--------------------------------------
113
## custom.covRob.Mcd <- function(R, ...){
114
## 
115
##   if(hasArg(control)) control = match.call(expand.dots = TRUE)$control else control = MycovRobMcd()
116
##   if(hasArg(alpha)) alpha = match.call(expand.dots = TRUE)$alpha else alpha = control$alpha
117
##   if(hasArg(nsamp)) nsamp = match.call(expand.dots = TRUE)$nsamp else nsamp = control$nsamp
118
##   if(hasArg(nmini)) nmini = match.call(expand.dots = TRUE)$nmini else nmini = control$nmini
119
##   if(hasArg(kmini)) kmini = match.call(expand.dots = TRUE)$kmini else kmini = control$kmini
120
##   if(hasArg(scalefn)) scalefn = match.call(expand.dots = TRUE)$scalefn else scalefn = control$scalefn
121
##   if(hasArg(maxcsteps)) maxcsteps = match.call(expand.dots = TRUE)$maxcsteps
122
##   else maxcsteps = control$maxcsteps
123
## 
124
##   if(hasArg(initHsets)) initHsets = match.call(expand.dots = TRUE)$initHsets
125
##   else initHsets = control$initHsets
126
## 
127
##   if(hasArg(seed)) seed = match.call(expand.dots = TRUE)$seed else seed = control$seed
128
##   if(hasArg(tolSolve)) tolSolve = match.call(expand.dots = TRUE)$tolSolve else tolSolve = control$tolSolve
129
##   if(hasArg(wgtFUN)) wgtFUN = match.call(expand.dots = TRUE)$wgtFUN else wgtFUN = control$wgtFUN
130
## 
131
##   if(hasArg(use.correction)) use.correction = match.call(expand.dots = TRUE)$use.correction
132
##   else use.correction = control$use.correction
133
## 
134
##   robustMCD <- robustbase::covMcd(x = R, alpha = alpha,
135
##                                   nsamp = nsamp, nmini = nmini,
136
##                                   kmini = kmini, seed = seed,
137
##                                   tolSolve = tolSolve, scalefn = scalefn,
138
##                                   maxcsteps = maxcsteps,
139
##                                   initHsets = initHsets,
140
##                                   wgtFUN = wgtFUN, use.correction = use.correction)
141
## 
142
##   return(list(mu = robustMCD$center, sigma = robustMCD$cov))
143
## }
144

145

146
## -------------------------------------------------------
147
opt <- optimize.portfolio(returns10, pspec, 
148
                   optimize_method = "CVXR", 
149
                   momentFUN = "custom.covRob.Mcd", 
150
                   alpha = 0.75, nsamp = 600)
151
outCovRobMcd <- opt.outputMvo(opt, returns10, digits = 3, frequency = "weekly")
152
(WtsCovRobMcd <- outCovRobMcd$Wgts)
153

154

155
## ---- eval = FALSE--------------------------------------
156
## MycovRobMcd <- function(alpha = 1/2,
157
##            nsamp = 500, nmini = 300, kmini = 5,
158
##            scalefn = "hrv2012", maxcsteps = 200,
159
##            seed = NULL, tolSolve = 1e-14,
160
##            wgtFUN = "01.original", beta,
161
##            use.correction = TRUE
162
##   ){
163
##     if(missing(beta) || !is.numeric(beta))
164
##       beta <- 0.975
165
## 
166
##     return(list(alpha = alpha, nsamp = nsamp,
167
##                 nmini = as.integer(nmini), kmini = as.integer(kmini),
168
##                 seed = as.integer(seed),
169
##                 tolSolve = tolSolve, scalefn = scalefn,
170
##                 maxcsteps = as.integer(maxcsteps),
171
##                 wgtFUN = wgtFUN, beta = beta,
172
##                 use.correction = use.correction))
173
##   }
174

175

176
## ---- eval = FALSE--------------------------------------
177
## covMcd.params <- MycovRobMcd(alpha = 0.75, nsamp = 600)
178
## opt <- optimize.portfolio(returns, pspec,
179
##                    optimize_method = "CVXR",
180
##                    momentFUN = "custom.covRob.Mcd",
181
##                    control = covMcd.params)
182

183

184
## ---- eval = FALSE--------------------------------------
185
## custom.covRob.TSGS <- function(R, ...){
186
##  if(hasArg(control)) control = match.call(expand.dots = TRUE)$control else control = MycovRobTSGS()
187
##  if(hasArg(filter)) filter = match.call(expand.dots = TRUE)$filter else filter = control$filter
188
## 
189
##  if(hasArg(partial.impute)) partial.impute = match.call(expand.dots = TRUE)$partial.impute
190
##  else partial.impute = control$partial.impute
191
## 
192
##  if(hasArg(tol)) tol = match.call(expand.dots=TRUE)$tol else tol=control$tol
193
##  if(hasArg(maxiter)) maxiter = match.call(expand.dots = TRUE)$maxiter else maxiter = control$maxiter
194
##  if(hasArg(loss)) loss = match.call(expand.dots = TRUE)$loss else loss = control$loss
195
##  if(hasArg(init)) init = match.call(expand.dots = TRUE)$init else init = control$init
196
## 
197
##   tsgsRob <- GSE::TSGS(x = R, filter = filter,
198
##                        partial.impute = partial.impute, tol = tol,
199
##                        maxiter = maxiter, method = loss,
200
##                        init = init)
201
## 
202
##   return(list(mu = tsgsRob@mu, sigma = tsgsRob@S))
203
## 
204
## }
205

206

207
## -------------------------------------------------------
208
opt <- optimize.portfolio(returns10, pspec, 
209
                   optimize_method = "CVXR", 
210
                   momentFUN = "custom.covRob.TSGS")
211
outCovRobTSGS <- opt.outputMvo(opt, returns10, digits = 3, frequency = "weekly")
212
(WtsCovRobTSGS <- outCovRobTSGS$Wgts)
213

214

215
## ---- eval = F------------------------------------------
216
## MycovRobTSGS <- function(filter = c("UBF-DDC","UBF","DDC","UF"),
217
##                 partial.impute = FALSE, tol = 1e-4, maxiter = 150,
218
##                 loss = c("bisquare","rocke"),
219
##                 init = c("emve", "qc", "huber", "imputed", "emve_c")){
220
## 
221
##   filter <- match.arg(filter)
222
##   loss <- match.arg(loss)
223
##   init <- match.arg(init)
224
## 
225
##   return(list(filter = filter, partial.impute = partial.impute,
226
##               tol = tol, maxiter = as.integer(maxiter),
227
##               loss = loss,init))
228
## }
229

230

231
## -------------------------------------------------------
232
dat <- data.frame(rbind(WtsCovClassic, WtsCovRobMM, WtsCovRobRocke,
233
                 WtsCovRobMcd, WtsCovRobTSGS))
234
print.data.frame(dat)
235

236

237
## -------------------------------------------------------
238
dat.mat <- rbind(outCovClassic[2:4], outCovRobMM[2:4],  outCovRobRocke[2:4],
239
                 outCovRobMcd[2:4], outCovRobTSGS[2:4])
240
dat <- as.data.frame(dat.mat)
241
row.names(dat) <- c("GmvLOcovClassic", "GmvLOcovRobMM", "GmvLOcovRobRocke",
242
                    "GmvLOcovRobMcd", "GmvLOcovRobTSGS")
243
print.data.frame(dat)
244

245

246
## -------------------------------------------------------
247
# Plot function
248
robPlot <- function(GMV, MARKET, plot=1){
249
  # Optimize Portfolio at Monthly Rebalancing and 5-Year Training
250
  if(plot == 1){
251
    momentEstFun = 'custom.covRob.Rocke'
252
    name = "GmvLOCovRobRocke"
253
  }else if(plot == 2){
254
    momentEstFun = 'custom.covRob.Mcd'
255
    name = "GmvLOCovMcd"
256
  }else if(plot == 3){
257
    momentEstFun = 'custom.covRob.TSGS'
258
    name = "GmvLOTSGS"
259
  }else{
260
    print("plot should be 1, 2 or 3")
261
    return()
262
  }
263
  
264
  bt.gmv.rob <- optimize.portfolio.rebalancing(returns, pspec,
265
                                            optimize_method = "CVXR",
266
                                            rebalance_on = "weeks",
267
                                            training_period = 100, 
268
                                            momentFUN = momentEstFun)
269
  
270
  # Extract time series of portfolio weights
271
  wts.gmv.rob <- extractWeights(bt.gmv.rob)
272
  # Compute cumulative returns of portfolio
273
  GMV.rob <- Return.rebalancing(returns, wts.gmv.rob)
274
  
275
  # Combine GMV.LO and MARKET cumulative returns
276
  ret.comb <- na.omit(merge(GMV.rob, GMV, MARKET, all=F))
277
  names(ret.comb) <- c(name, "GmvLO", "MARKET")
278
  
279
  plot <- backtest.plot(ret.comb, colorSet = c("darkgreen", "black", 
280
                        "red"), ltySet = c(1,2,3))
281
  
282
  return(list(ret = ret.comb, plot = plot))
283
}
284

285

286

287
## -------------------------------------------------------
288
bt.gmv <- optimize.portfolio.rebalancing(returns, pspec,
289
                                         optimize_method = "CVXR",
290
                                         rebalance_on="weeks",
291
                                         training_period = 100)
292
# Extract time series of portfolio weights
293
wts.gmv <- extractWeights(bt.gmv)
294
# Compute cumulative returns of portfolio
295
GMV <- Return.rebalancing(returns, wts.gmv)
296
  
297
res.covRob <- robPlot(GMV=GMV, MARKET=MARKET, plot=1)
298
res.covRob$plot
299

300

301
## -------------------------------------------------------
302
table.Drawdowns(res.covRob$ret$GmvLOCovRobRocke, top=1)
303

304

305
## -------------------------------------------------------
306
table.Drawdowns(res.covRob$ret$GmvLO, top=1)
307

308

309
## -------------------------------------------------------
310
table.Drawdowns(res.covRob$ret$MARKET, top=1)
311

312

313
## -------------------------------------------------------
314
set.seed(1234)
315
res.covMcd = robPlot(GMV=GMV, MARKET=MARKET, plot=2)
316
res.covMcd$plot
317

318

319
## -------------------------------------------------------
320
# longest drawdown for robust based portfolio
321
table.Drawdowns(res.covMcd$ret$GmvLOCovMcd, top=1) 
322

323

324
## -------------------------------------------------------
325
res.TSGS = robPlot(GMV=GMV, MARKET=MARKET, plot=3)
326
res.TSGS$plot
327

328

329
## -------------------------------------------------------
330
# longest drawdown for robust based portfolio
331
table.Drawdowns(res.TSGS$ret$GmvLOTSGS, top=1)
332

333

334
## ---- eval = FALSE--------------------------------------
335
## user.covRob.Rocke <- function(R){
336
##   out <- list()
337
##   robustCov <- RobStatTM::covRobRocke(X = R, maxit = 200, maxsteps = 10)
338
##   out$sigma <- robustCov$cov
339
##   out$mu <- robustCov$center
340
##   return(out)
341
## }
342

343

344
## ---- eval = FALSE--------------------------------------
345
## opt <- optimize.portfolio(returns, pspec,
346
##                    optimize_method = "CVXR",
347
##                    momentFUN = "user.covRob.Rocke")
348

349

350
## ---- eval = FALSE--------------------------------------
351
## user.covRob.OGK <- function(R){
352
##   robustCov <- robustbase::covOGK(x = R)
353
##   return(list(mu = robustCov$center, sigma = robustCov$cov))
354
## }
355

356

357
## ---- eval = FALSE--------------------------------------
358
## opt <- optimize.portfolio(returns, pspec,
359
##                    optimize_method = "CVXR",
360
##                    momentFUN = "user.covMcd.OGK")
361

362
## opt$moment_values
363

364
Product

Resources

Company
