-
Notifications
You must be signed in to change notification settings - Fork 16
/
Copy pathFRE7241_Lecture_2.R
1910 lines (1910 loc) · 77.2 KB
/
FRE7241_Lecture_2.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# Create a design matrix of IEF and VTI returns
desm <- na.omit(rutils::etfenv$returns[, c("IEF", "VTI")])
retvti <- desm$VTI
# Add returns with perfect timing skill
desm <- cbind(desm, 0.5*(retvti+abs(retvti)), retvti^2)
colnames(desm)[3:4] <- c("merton", "treynor")
# Perform Merton-Henriksson test regression
regmod <- lm(IEF ~ VTI + merton, data=desm); summary(regmod)
# Perform Treynor-Mazuy test regression
regmod <- lm(IEF ~ VTI + treynor, data=desm); summary(regmod)
# Plot residual scatterplot
resids <- (desm$IEF - regmod$coeff["VTI"]*retvti)
plot.default(x=retvti, y=resids, xlab="VTI", ylab="IEF")
title(main="Treynor-Mazuy Market Timing Test\n for IEF vs VTI", line=0.5)
# Plot fitted (predicted) response values
coefreg <- summary(regmod)$coeff
fitv <- regmod$fitted.values - coefreg["VTI", "Estimate"]*retvti
tvalue <- round(coefreg["treynor", "t value"], 2)
points.default(x=retvti, y=fitv, pch=16, col="red")
text(x=0.0, y=0.8*max(resids), paste("Treynor test t-value =", tvalue))
library(rutils)
# Extract the ETF prices from rutils::etfenv$prices
pricev <- rutils::etfenv$prices
pricev <- zoo::na.locf(pricev, na.rm=FALSE)
pricev <- zoo::na.locf(pricev, fromLast=TRUE)
datev <- zoo::index(pricev)
# Calculate the dollar returns
retd <- rutils::diffit(pricev)
# Or
# retd <- lapply(pricev, rutils::diffit)
# retd <- rutils::do_call(cbind, retd)
# Calculate the percentage returns
retp <- retd/rutils::lagit(pricev, lagg=1, pad_zeros=FALSE)
# Calculate the log returns
retl <- rutils::diffit(log(pricev))
# Set the initial dollar returns
retd[1, ] <- pricev[1, ]
# Calculate the prices from dollar returns
pricen <- cumsum(retd)
all.equal(pricen, pricev)
# Compound the percentage returns
pricen <- cumprod(1 + retp)
# Set the initial prices
prici <- as.numeric(pricev[1, ])
pricen <- lapply(1:NCOL(pricen), function (i) prici[i]*pricen[, i])
pricen <- rutils::do_call(cbind, pricen)
# pricen <- t(t(pricen)*prici)
all.equal(pricen, pricev, check.attributes=FALSE)
# Plot log VTI prices
endd <- rutils::calc_endpoints(rutils::etfenv$VTI, interval="weeks")
dygraphs::dygraph(log(quantmod::Cl(rutils::etfenv$VTI)[endd]),
main="Logarithm of VTI Prices") %>%
dyOptions(colors="blue", strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Calculate the percentage VTI returns
pricev <- rutils::etfenv$prices$VTI
pricev <- na.omit(pricev)
retp <- rutils::diffit(pricev)/rutils::lagit(pricev, lagg=1, pad_zeros=FALSE)
# Funding rate per day
ratef <- 0.01/252
# Margin account value
marginv <- cumsum(retp)
# Cumulative funding costs
costf <- cumsum(ratef*marginv)
# Add funding costs to margin account
marginv <- (marginv + costf)
# dygraph plot of margin and funding costs
datav <- cbind(marginv, costf)
colnamev <- c("Margin", "Cumulative Funding")
colnames(datav) <- colnamev
endd <- rutils::calc_endpoints(datav, interval="weeks")
dygraphs::dygraph(datav[endd], main="VTI Margin Funding Costs") %>%
dyAxis("y", label=colnamev[1], independentTicks=TRUE) %>%
dyAxis("y2", label=colnamev[2], independentTicks=TRUE) %>%
dySeries(name=colnamev[1], axis="y", col="blue") %>%
dySeries(name=colnamev[2], axis="y2", col="red", strokeWidth=3) %>%
dyLegend(show="always", width=300)
# bidask equal to 1 bp for liquid ETFs
bidask <- 0.001
# Cumulative transaction costs
costs <- bidask*cumsum(abs(retp))/2
# Subtract transaction costs from margin account
marginv <- cumsum(retp)
marginv <- (marginv - costs)
# dygraph plot of margin and transaction costs
datav <- cbind(marginv, costs)
colnamev <- c("Margin", "Cumulative Transaction Costs")
colnames(datav) <- colnamev
dygraphs::dygraph(datav[endd], main="VTI Transaction Costs") %>%
dyAxis("y", label=colnamev[1], independentTicks=TRUE) %>%
dyAxis("y2", label=colnamev[2], independentTicks=TRUE) %>%
dySeries(name=colnamev[1], axis="y", col="blue") %>%
dySeries(name=colnamev[2], axis="y2", col="red", strokeWidth=3) %>%
dyLegend(show="always", width=300)
# Calculate the VTI and IEF dollar returns
pricev <- rutils::etfenv$prices[, c("VTI", "IEF")]
pricev <- na.omit(pricev)
retd <- rutils::diffit(pricev)
datev <- zoo::index(pricev)
# Calculate the VTI and IEF percentage returns
retp <- retd/rutils::lagit(pricev, lagg=1, pad_zeros=FALSE)
# Wealth of fixed shares equal to $0.5 each at start (without rebalancing)
weightv <- c(0.5, 0.5) # dollar weights
wealthfs <- drop(cumprod(1 + retp) %*% weightv)
# Or using the dollar returns
prici <- as.numeric(pricev[1, ])
retd[1, ] <- pricev[1, ]
wealthfs2 <- cumsum(retd %*% (weightv/prici))
all.equal(wealthfs, drop(wealthfs2))
# Wealth of fixed dollars equal to $0.5 each (with rebalancing)
wealthfd <- cumsum(retp %*% weightv)
# Calculate the Sharpe and Sortino ratios
wealthv <- cbind(log(wealthfs), wealthfd)
wealthv <- xts::xts(wealthv, datev)
colnames(wealthv) <- c("Fixed shares", "Fixed dollars")
sqrt(252)*sapply(rutils::diffit(wealthv), function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# Plot the log wealth
colnamev <- colnames(wealthv)
endd <- rutils::calc_endpoints(retp, interval="weeks")
dygraphs::dygraph(wealthv[endd], main="Wealth of Weighted Portfolios") %>%
dySeries(name=colnamev[1], col="blue", strokeWidth=2) %>%
dySeries(name=colnamev[2], col="red", strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Margin account for fixed dollars (with rebalancing)
marginv <- cumsum(retp %*% weightv)
# Cumulative transaction costs
costs <- bidask*cumsum(abs(retp) %*% weightv)/2
# Subtract transaction costs from margin account
marginv <- (marginv - costs)
# dygraph plot of margin and transaction costs
datav <- cbind(marginv, costs)
datav <- xts::xts(datav, datev)
colnamev <- c("Margin", "Cumulative Transaction Costs")
colnames(datav) <- colnamev
dygraphs::dygraph(datav[endd], main="Fixed Dollar Portfolio Transaction Costs") %>%
dyAxis("y", label=colnamev[1], independentTicks=TRUE) %>%
dyAxis("y2", label=colnamev[2], independentTicks=TRUE) %>%
dySeries(name=colnamev[1], axis="y", col="blue") %>%
dySeries(name=colnamev[2], axis="y2", col="red", strokeWidth=3) %>%
dyLegend(show="always", width=300)
# Wealth of fixed shares (without rebalancing)
wealthfs <- cumsum(retd %*% (weightv/prici))
# Or compound the percentage returns
wealthfs <- drop(cumprod(1 + retp) %*% weightv)
# Wealth of equal wealth strategy (with rebalancing)
wealthew <- cumprod(1 + retp %*% weightv)
wealthv <- cbind(wealthfs, wealthew)
wealthv <- xts::xts(wealthv, datev)
colnames(wealthv) <- c("Fixed shares", "Prop dollars")
wealthv <- log(wealthv)
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(rutils::diffit(wealthv), function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# Plot the log wealth
dygraphs::dygraph(wealthv[endd],
main="Wealth of Proportional Wealth Allocations") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Returns in excess of weighted returns
retw <- retp %*% weightv
retx <- lapply(retp, function(x) (retw - x))
retx <- do.call(cbind, retx)
sum(retx %*% weightv)
# Calculate the weighted sum of absolute excess returns
retx <- abs(retx) %*% weightv
# Total dollar amount of stocks that need to be traded
retx <- retx*rutils::lagit(wealthew)
# Cumulative transaction costs
costs <- bidask*cumsum(retx)/2
# Subtract transaction costs from wealth
wealthew <- (wealthew - costs)
# dygraph plot of wealth and transaction costs
wealthv <- cbind(wealthew, costs)
wealthv <- xts::xts(wealthv, datev)
colnamev <- c("Wealth", "Cumulative Transaction Costs")
colnames(wealthv) <- colnamev
dygraphs::dygraph(wealthv[endd],
main="Transaction Costs With Equal Wealths") %>%
dyAxis("y", label=colnamev[1], independentTicks=TRUE) %>%
dyAxis("y2", label=colnamev[2], independentTicks=TRUE) %>%
dySeries(name=colnamev[1], axis="y", col="blue") %>%
dySeries(name=colnamev[2], axis="y2", col="red", strokeWidth=3) %>%
dyLegend(show="always", width=300)
# Wealth of fixed shares (without rebalancing)
wealthfs <- drop(apply(retp, 2, function(x) cumprod(1 + x)) %*% weightv)-1
# Wealth of proportional wealth allocations (with rebalancing)
wealthew <- cumprod(1 + retp %*% weightv) - 1
# Wealth of proportional target allocation (with rebalancing)
retp <- zoo::coredata(retp)
threshv <- 0.05
wealthv <- matrix(nrow=NROW(retp), ncol=2)
colnames(wealthv) <- colnames(retp)
wealthv[1, ] <- weightv
for (it in 2:NROW(retp)) {
# Accrue wealth without rebalancing
wealthv[it, ] <- wealthv[it-1, ]*(1 + retp[it, ])
# Rebalance if wealth allocations differ from weights
if (sum(abs(wealthv[it, ] - sum(wealthv[it, ])*weightv))/sum(wealthv[it, ]) > threshv) {
# cat("Rebalance at:", it, "\n")
wealthv[it, ] <- sum(wealthv[it, ])*weightv
} # end if
} # end for
wealthv <- rowSums(wealthv) - 1
wealthv <- cbind(wealthew, wealthv)
wealthv <- xts::xts(wealthv, datev)
colnames(wealthv) <- c("Equal Wealths", "Proportional Target")
dygraphs::dygraph(wealthv, main="Wealth of Proportional Target Allocations") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Calculate the stock and bond returns
retp <- na.omit(rutils::etfenv$returns[, c("VTI", "IEF")])
weightv <- c(0.4, 0.6)
retp <- cbind(retp, retp %*% weightv)
colnames(retp)[3] <- "Combined"
# Calculate the correlations
cor(retp)
# Calculate the Sharpe ratios
sqrt(252)*sapply(retp, function(x) mean(x)/sd(x))
# Calculate the standard deviation, skewness, and kurtosis
sapply(retp, function(x) {
# Calculate the standard deviation
stdev <- sd(x)
# Standardize the returns
x <- (x - mean(x))/stdev
c(stdev=stdev, skew=mean(x^3), kurt=mean(x^4))
}) # end sapply
# Wealth of equal wealth strategy
wealthv <- cumsum(retp)
# Calculate the a vector of monthly end points
endd <- rutils::calc_endpoints(retp, interval="weeks")
# Plot cumulative log wealth
dygraphs::dygraph(wealthv[endd],
main="Stocks and Bonds With Equal Wealths") %>%
dyOptions(colors=c("blue", "green", "blue", "red")) %>%
dySeries("Combined", color="red", strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Calculate the Sharpe ratios
sqrt(252)*sapply(retp, function(x) mean(x)/sd(x))
# Calculate the Sharpe ratios for vector of weights
weightv <- seq(0.05, 0.95, 0.05)
sharpev <- sqrt(252)*sapply(weightv, function(weight) {
weightv <- c(weight, 1-weight)
retp <- (retp[, 1:2] %*% weightv)
mean(retp)/sd(retp)
}) # end sapply
# Calculate the optimal VTI weight
weightm <- weightv[which.max(sharpev)]
# Calculate the optimal weight using optimization
calc_sharpe <- function(weight) {
weightv <- c(weight, 1-weight)
retp <- (retp[, 1:2] %*% weightv)
-mean(retp)/sd(retp)
} # end calc_sharpe
optv <- optimize(calc_sharpe, interval=c(0, 1))
weightm <- optv$minimum
# Plot Sharpe ratios
plot(x=weightv, y=sharpev,
main="Sharpe Ratio as Function of VTI Weight",
xlab="VTI weight", ylab="Sharpe Ratio",
t="l", lwd=3, col="blue")
abline(v=weightm, lty="dashed", lwd=1, col="blue")
text(x=weightm, y=0.7*max(sharpev), pos=4, cex=1.2,
labels=paste("optimal VTI weight =", round(weightm, 2)))
# Coerce the returns from xts time series to matrix
retp <- zoo::coredata(retp[, 1:2])
nrows <- NROW(retp)
# Bootstrap the returns and Calculate the a list of random returns
nboot <- 1e4
library(parallel) # Load package parallel
ncores <- detectCores() - 1 # Number of cores
# Perform parallel bootstrap under Windows
compclust <- makeCluster(ncores) # Initialize compute cluster under Windows
clusterSetRNGStream(compclust, 1121) # Reset random number generator in all cores
clusterExport(compclust, c("retp", "nrows"))
bootd <- parLapply(compclust, 1:nboot, function(x) {
retp[sample.int(nrows, replace=TRUE), ]
}) # end parLapply
# Perform parallel bootstrap under Mac-OSX or Linux
set.seed(1121, "Mersenne-Twister", sample.kind="Rejection")
bootd <- mclapply(1:nboot, function(x) {
retp[sample.int(nrows, replace=TRUE), ]
}, mc.cores=ncores) # end mclapply
is.list(bootd); NROW(bootd); dim(bootd[[1]])
# Calculate the distribution of terminal wealths under Windows
wealthv <- parLapply(compclust, bootd, function(retp) {
apply(retp, 2, function(x) prod(1 + x))
}) # end parLapply
# Calculate the distribution of terminal wealths under Mac-OSX or Linux
wealthv <- mclapply(bootd, function(retp) {
apply(retp, 2, function(x) prod(1 + x))
}, mc.cores=ncores) # end mclapply
wealthv <- do.call(rbind, wealthv)
class(wealthv); dim(wealthv); tail(wealthv)
# Calculate the means and standard deviations of the terminal wealths
apply(wealthv, 2, mean)
apply(wealthv, 2, sd)
# Extract the terminal wealths of VTI and IEF
vtiw <- wealthv[, "VTI"]
iefw <- wealthv[, "IEF"]
# Plot the densities of the terminal wealths of VTI and IEF
vtim <- mean(vtiw); iefm <- mean(iefw)
vtid <- density(vtiw); iefd <- density(iefw)
plot(vtid, col="blue", lwd=3, xlab="wealth",
xlim=c(0, 2*max(iefd$x)), ylim=c(0, max(iefd$y)),
main="Terminal Wealth Distributions of VTI and IEF")
lines(iefd, col="green", lwd=3)
abline(v=vtim, col="blue", lwd=2, lty="dashed")
text(x=vtim, y=0.5, labels="VTI mean", pos=4, cex=0.8)
abline(v=iefm, col="green", lwd=2, lty="dashed")
text(x=iefm, y=0.5, labels="IEF mean", pos=4, cex=0.8)
legend(x="topright", legend=c("VTI", "IEF"),
inset=0.1, cex=1.0, bg="white", bty="n", y.intersp=0.5,
lwd=6, lty=1, col=c("blue", "green"))
# Calculate the distributions of stock wealth
holdv <- nrows*seq(0.1, 1.0, 0.1)
wealthm <- mclapply(bootd, function(retp) {
sapply(holdv, function(holdp) {
prod(1 + retp[1:holdp, "VTI"])
}) # end sapply
}, mc.cores=ncores) # end mclapply
wealthm <- do.call(rbind, wealthm)
dim(wealthm)
# Plot the stock wealth for long and short holding periods
wealth1 <- wealthm[, 9]
wealth2 <- wealthm[, 1]
mean1 <- mean(wealth1); mean2 <- mean(wealth2)
dens1 <- density(wealth1); dens2 <- density(wealth2)
plot(dens1, col="blue", lwd=3, xlab="wealth",
xlim=c(0, 3*max(dens2$x)), ylim=c(0, max(dens2$y)),
main="Wealth Distributions for Long and Short Holding Periods")
lines(dens2, col="green", lwd=3)
abline(v=mean1, col="blue", lwd=2, lty="dashed")
text(x=mean1, y=0.5, labels="Long", pos=4, cex=0.8)
abline(v=mean2, col="green", lwd=2, lty="dashed")
text(x=mean2, y=0.5, labels="Short", pos=4, cex=0.8)
legend(x="top", legend=c("Long", "Short"),
inset=0.1, cex=1.0, bg="white", bty="n", y.intersp=0.5,
lwd=6, lty=1, col=c("blue", "green"))
# Define the risk-adjusted wealth measure
riskretfun <- function(wealthv) {
riskv <- 0.01 # Risk floor
if (min(wealthv) < 1) # Some wealth is below par
# Calculate the mean stock wealth below par
riskv <- mean((1-wealthv)[wealthv<1])
mean(wealthv)/riskv
} # end riskretfun
# Calculate the stock wealth risk-return ratios
riskrets <- apply(wealthm, 2, riskretfun)
# Plot the stock wealth risk-return ratios
plot(x=holdv, y=riskrets,
main="Stock Risk-Return Ratio as Function of Holding Period",
xlab="Holding Period", ylab="Ratio",
t="l", lwd=3, col="blue")
# Calculate the distributions of portfolio wealth
weightv <- seq(0.05, 0.95, 0.05)
wealthm <- mclapply(bootd, function(retp) {
sapply(weightv, function(weight) {
prod(1 + retp %*% c(weight, 1-weight))
}) # end sapply
}, mc.cores=ncores) # end mclapply
wealthm <- do.call(rbind, wealthm)
dim(wealthm)
# Calculate the portfolio risk-return ratios
riskrets <- apply(wealthm, 2, riskretfun)
# Calculate the optimal VTI weight
weightm <- weightv[which.max(riskrets)]
# Plot the portfolio risk-return ratios
plot(x=weightv, y=riskrets,
main="Portfolio Risk-Return Ratio as Function of VTI Weight",
xlab="VTI weight", ylab="Ratio",
t="l", lwd=3, col="blue")
abline(v=weightm, lty="dashed", lwd=1, col="blue")
text(x=weightm, y=0.7*max(riskrets), pos=4, cex=1.2,
labels=paste("optimal VTI weight =", round(weightm, 2)))
# Extract the ETF returns
symbolv <- c("VTI", "IEF", "DBC")
retp <- na.omit(rutils::etfenv$returns[, symbolv])
# Calculate the all-weather portfolio wealth
weightaw <- c(0.30, 0.55, 0.15)
retp <- cbind(retp, retp %*% weightaw)
colnames(retp)[4] <- "All Weather"
# Calculate the Sharpe ratios
sqrt(252)*sapply(retp, function(x) mean(x)/sd(x))
# Calculate the cumulative wealth from returns
wealthv <- cumsum(retp)
# Calculate the a vector of monthly end points
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
# dygraph all-weather wealth
dygraphs::dygraph(wealthv[endd], main="All-Weather Portfolio") %>%
dyOptions(colors=c("blue", "green", "orange", "red")) %>%
dySeries("All Weather", color="red", strokeWidth=2) %>%
dyLegend(show="always", width=400)
# Plot all-weather wealth
plot_theme <- chart_theme()
plot_theme$col$line.col <- c("orange", "blue", "green", "red")
quantmod::chart_Series(wealthv, theme=plot_theme, lwd=c(2, 2, 2, 4),
name="All-Weather Portfolio")
legend("topleft", legend=colnames(wealthv),
inset=0.1, bg="white", lty=1, lwd=6, y.intersp=0.5,
col=plot_theme$col$line.col, bty="n")
# Calculate the VTI returns
retp <- na.omit(rutils::etfenv$returns$VTI["2008/2009"])
datev <- zoo::index(retp)
nrows <- NROW(retp)
retp <- drop(zoo::coredata(retp))
# Bond floor
bfloor <- 60
# CPPI multiplier
coeff <- 2
# Portfolio market values
portfv <- numeric(nrows)
# Initial principal
portfv[1] <- 100
# Stock allocation
stockv <- numeric(nrows)
stockv[1] <- min(coeff*(portfv[1] - bfloor), portfv[1])
# Bond allocation
bondv <- numeric(nrows)
bondv[1] <- (portfv[1] - stockv[1])
# Simulate CPPI strategy
for (t in 2:nrows) {
portfv[t] <- portfv[t-1] + stockv[t-1]*retp[t]
stockv[t] <- min(coeff*(portfv[t] - bfloor), portfv[t])
bondv[t] <- (portfv[t] - stockv[t])
} # end for
# dygraph plot of CPPI strategy
pricev <- 100*cumprod(1 + retp)
datav <- xts::xts(cbind(stockv, bondv, portfv, pricev), datev)
colnames(datav) <- c("stocks", "bonds", "CPPI", "VTI")
endd <- rutils::calc_endpoints(datav, interval="weeks")
dygraphs::dygraph(datav[endd], main="CPPI strategy") %>%
dyOptions(colors=c("red", "green", "blue", "orange"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Calculate the dollar returns of VTI and IEF
pricev <- na.omit(rutils::etfenv$prices[, c("VTI", "IEF")])
retd <- rutils::diffit(pricev)
# Scale the stock prices to $1 at beginning
prici <- as.numeric(pricev[1, ]) # Initial stock prices
pricesc <- pricev
pricesc$VTI <- pricesc$VTI/prici[1]
pricesc$IEF <- pricesc$IEF/prici[2]
sum(pricesc[1, ])
retsc <- rutils::diffit(pricesc)
# Wealth of fixed number of shares (without rebalancing)
weightv <- c(0.5, 0.5) # Buy $0.5 of each stock
wealthed <- 1 + cumsum(retsc %*% weightv)
# Calculate the stock prices with unit dollar volatility
stdev <- sapply(retd, sd)
pricesd <- pricev
pricesd$VTI <- pricev$VTI/stdev["VTI"]
pricesd$IEF <- pricev$IEF/stdev["IEF"]
retsd <- rutils::diffit(pricesd)
sapply(retsd, sd)
# Or equivalent using lapply()
# Calculate the standardized dollar returns
retsd <- lapply(retd, function(x) x/sd(x))
retsd <- do.call(cbind, retsd)
sapply(retsd, sd)
# Initial stock prices with unit dollar volatility
prici <- as.numeric(pricev[1, ])/sapply(retd, sd)
# Wealth of shares with equal dollar volatilities
wealthev <- 1 + cumsum(retsd %*% weightv)/sum(prici*weightv)
# Scale the sum of stock prices at beginning to $1
pricesd <- pricesd/sum(pricesd[1, ])
retsd <- rutils::diffit(pricesd)
sapply(retsd, sd)
# Wealth of shares with equal dollar volatilities
wealthev <- 1 + cumsum(retsd %*% weightv)
# Calculate the Sharpe and Sortino ratios
wealthv <- xts::xts(cbind(wealthed, wealthev), zoo::index(pricev))
colnames(wealthv) <- c("Equal dollar", "Equal volatility")
sqrt(252)*sapply(rutils::diffit(wealthv), function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# Plot the log wealth
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(log(wealthv[endd]),
main="Wealth of Equal Dollar And Equal Volatility") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Calculate the wealth of proportional wealth allocations (with rebalancing)
retp <- retd/rutils::lagit(pricev, lagg=1, pad_zeros=FALSE)
weightv <- c(0.5, 0.5)
wealthew <- cumprod(1 + retp %*% weightv)
# Calculate the trailing dollar volatilities
vold <- HighFreq::run_var(retd, lambda=0.9)
vold <- sqrt(vold)
vold[vold == 0] <- 1
# Calculate the standardized prices with unit dollar volatilities
pricerp <- pricev/vold
pricerp <- rutils::lagit(pricerp)
pricerp[1, ] <- 1
# Scale the sum of stock prices to $1
pricerp <- pricerp/rowSums(pricerp)
# Calculate the percentage returns of risk parity
retrp <- retp*pricerp
# Calculate the wealth of risk parity
wealthrp <- cumprod(1 + retrp %*% weightv)
# Calculate the log wealths
wealthv <- cbind(wealthew, wealthrp)
wealthv <- xts::xts(wealthv, zoo::index(pricev))
colnames(wealthv) <- c("PropDollars", "Risk Parity")
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(rutils::diffit(wealthv), function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# Plot a dygraph of the log wealths
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(log(wealthv[endd]),
main="Log Wealth of Risk Parity vs Equal Wealths") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Test risk parity market timing of VTI using Treynor-Mazuy test
retrp <- rutils::diffit(wealthv)
retvti <- retp$VTI
desm <- cbind(retrp, retvti, retvti^2)
colnames(desm)[1:2] <- c("prop", "riskp")
colnames(desm)[4] <- "Treynor"
regmod <- lm(riskp ~ VTI + Treynor, data=desm)
summary(regmod)
# Plot residual scatterplot
resids <- regmod$residuals
plot.default(x=retvti, y=resids, xlab="VTI", ylab="residuals")
title(main="Treynor-Mazuy Market Timing Test\n for Risk Parity vs VTI", line=0.5)
# Plot fitted (predicted) response values
coefreg <- summary(regmod)$coeff
fitv <- regmod$fitted.values - coefreg["VTI", "Estimate"]*retvti
tvalue <- round(coefreg["Treynor", "t value"], 2)
points.default(x=retvti, y=fitv, pch=16, col="red")
text(x=0.0, y=0.8*max(resids), paste("Treynor test t-value =", tvalue))
# Test for equal wealth strategy market timing of VTI using Treynor-Mazuy test
regmod <- lm(prop ~ VTI + Treynor, data=desm)
summary(regmod)
# Plot fitted (predicted) response values
coefreg <- summary(regmod)$coeff
fitv <- regmod$fitted.values - coefreg["VTI", "Estimate"]*retvti
points.default(x=retvti, y=fitv, pch=16, col="blue")
text(x=0.0, y=0.6*max(resids), paste("Prop Alloc t-value =", round(coefreg["Treynor", "t value"], 2)))
# Calculate the positions
retp <- na.omit(rutils::etfenv$returns$VTI)
posv <- rep(NA_integer_, NROW(retp))
datev <- zoo::index(retp)
datev <- format(datev, "%m-%d")
posv[datev == "05-01"] <- 0
posv[datev == "05-03"] <- 0
posv[datev == "11-01"] <- 1
posv[datev == "11-03"] <- 1
# Carry forward and backward non-NA posv
posv <- zoo::na.locf(posv, na.rm=FALSE)
posv <- zoo::na.locf(posv, fromLast=TRUE)
# Calculate the strategy returns
pnlinmay <- posv*retp
wealthv <- cbind(retp, pnlinmay)
colnames(wealthv) <- c("VTI", "sell_in_may")
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(wealthv, function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# Plot wealth of Sell in May strategy
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(cumsum(wealthv)[endd], main="Sell in May Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# OR: Open x11 for plotting
x11(width=6, height=5)
par(mar=c(4, 4, 3, 1), oma=c(0, 0, 0, 0))
plot_theme <- chart_theme()
plot_theme$col$line.col <- c("blue", "red")
quantmod::chart_Series(wealthv, theme=plot_theme, name="Sell in May Strategy")
legend("topleft", legend=colnames(wealthv),
inset=0.1, bg="white", lty=1, lwd=6, y.intersp=0.5,
col=plot_theme$col$line.col, bty="n")
# Test if Sell in May strategy can time VTI
desm <- cbind(wealthv, 0.5*(retp+abs(retp)), retp^2)
colnames(desm)[3:4] <- c("Merton", "Treynor")
# Perform Merton-Henriksson test
regmod <- lm(sell_in_may ~ VTI + Merton, data=desm)
summary(regmod)
# Perform Treynor-Mazuy test
regmod <- lm(sell_in_may ~ VTI + Treynor, data=desm)
summary(regmod)
# Plot Treynor-Mazuy residual scatterplot
resids <- (pnlinmay - regmod$coeff["VTI"]*retp)
plot.default(x=retp, y=resids, xlab="VTI", ylab="residuals")
title(main="Treynor-Mazuy Market Timing Test\n for Sell in May vs VTI", line=0.5)
# Plot fitted (predicted) response values
coefreg <- summary(regmod)$coeff
fitv <- regmod$fitted.values - coefreg["VTI", "Estimate"]*retp
tvalue <- round(coefreg["Treynor", "t value"], 2)
points.default(x=retp, y=fitv, pch=16, col="red")
text(x=0.0, y=0.8*max(resids), paste("Treynor test t-value =", tvalue))
# Calculate the log of OHLC VTI prices
ohlc <- log(rutils::etfenv$VTI)
openp <- quantmod::Op(ohlc)
highp <- quantmod::Hi(ohlc)
lowp <- quantmod::Lo(ohlc)
closep <- quantmod::Cl(ohlc)
# Calculate the close-to-close log returns,
# the daytime open-to-close returns
# and the overnight close-to-open returns.
retp <- rutils::diffit(closep)
colnames(retp) <- "daily"
retd <- (closep - openp)
colnames(retd) <- "daytime"
reton <- (openp - rutils::lagit(closep, lagg=1, pad_zeros=FALSE))
colnames(reton) <- "overnight"
# Calculate the Sharpe and Sortino ratios
wealthv <- cbind(retp, reton, retd)
sqrt(252)*sapply(wealthv, function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# Plot the log wealth
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(cumsum(wealthv)[endd],
main="Wealth of Close-to-Close, Overnight, and Daytime Strategies") %>%
dySeries(name="daily", strokeWidth=2, col="blue") %>%
dySeries(name="overnight", strokeWidth=2, col="red") %>%
dySeries(name="daytime", strokeWidth=2, col="green") %>%
dyLegend(width=500)
# Calculate the VTI returns
retp <- na.omit(rutils::etfenv$returns$VTI)
datev <- zoo::index(retp)
# Calculate the first business day of every month
dayv <- as.numeric(format(datev, "%d"))
indeks <- which(rutils::diffit(dayv) < 0)
datev[head(indeks)]
# Calculate the Turn of the Month dates
indeks <- lapply((-1):2, function(x) indeks + x)
indeks <- do.call(c, indeks)
sum(indeks > NROW(datev))
indeks <- sort(indeks)
datev[head(indeks, 11)]
# Calculate the Turn of the Month pnls
pnls <- numeric(NROW(retp))
pnls[indeks] <- retp[indeks, ]
# Combine data
wealthv <- cbind(retp, pnls)
colnamev <- c("VTI", "TOM Strategy")
colnames(wealthv) <- colnamev
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(wealthv, function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# dygraph plot VTI Turn of the Month strategy
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(cumsum(wealthv)[endd],
main="Turn of the Month Strategy") %>%
dyAxis("y", label=colnamev[1], independentTicks=TRUE) %>%
dyAxis("y2", label=colnamev[2], independentTicks=TRUE) %>%
dySeries(name=colnamev[1], axis="y", strokeWidth=2, col="blue") %>%
dySeries(name=colnamev[2], axis="y2", strokeWidth=2, col="red")
# Calculate the VTI prices and returns
pricev <- na.omit(rutils::etfenv$prices$VTI)
nrows <- NROW(pricev)
datev <- zoo::index(pricev)
retp <- rutils::diffit(log(pricev))
# Simulate stop-loss strategy
stopl <- 0.05 # Stop-loss percentage
pricem <- cummax(pricev) # Trailing maximum prices
# Calculate the drawdown
dd <- (pricev - pricem)
pnls <- retp # Initialize PnLs
for (i in 1:(nrows-1)) {
# Check for stop-loss
if (dd[i] < -stopl*pricem[i])
pnls[i+1] <- 0 # Set PnLs = 0 if in stop-loss
} # end for
# Same but without using loops in R
pnls2 <- retp
insl <- rutils::lagit(dd < -stopl*pricem)
pnls2 <- ifelse(insl, 0, pnls2)
all.equal(pnls, pnls2, check.attributes=FALSE)
# Combine the data
wealthv <- cbind(retp, pnls)
colnamev <- c("VTI", "Strategy")
colnames(wealthv) <- colnamev
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(wealthv, function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# dygraph plot the stop-loss strategy
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(cumsum(wealthv)[endd],
main="VTI Stop-loss Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Plot dygraph with shading
# Create colors for background shading
indic <- (rutils::diffit(insl) != 0) # Indices of stop-loss
crossd <- c(datev[indic], datev[nrows]) # Dates of stop-loss
shadev <- ifelse(insl[indic] == 1, "antiquewhite", "lightgreen")
# Create dygraph object without plotting it
dyplot <- dygraphs::dygraph(cumsum(wealthv),
main="VTI Stop-loss Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Add shading to dygraph object
for (i in 1:NROW(shadev)) {
dyplot <- dyplot %>% dyShading(from=crossd[i], to=crossd[i+1], color=shadev[i])
} # end for
# Plot the dygraph object
dyplot
# Simulate multiple stop-loss strategies
dd <- (pricev - pricem)
stopv <- 0.01*(1:30)
pnlc <- sapply(stopv, function(stopl) {
pnls <- retp
insl <- rutils::lagit(dd < -stopl*pricem)
pnls <- ifelse(insl, 0, pnls)
sum(pnls)
}) # end sapply
# Plot cumulative pnls for stop-loss strategies
plot(x=stopv, y=pnlc,
main="Cumulative PnLs for Stop-loss Strategies",
xlab="stop-loss percent", ylab="cumulative pnl",
t="l", lwd=3, col="blue")
# Define function for simulating a stop-start strategy
sim_stopstart <- function(stopl) {
maxp <- pricev[1] # Trailing maximum price
minp <- pricev[1] # Trailing minimum price
insl <- FALSE # Is in stop-loss?
insg <- FALSE # Is in start-gain?
pnls <- retp # Initialize PnLs
for (i in 1:nrows) {
if (insl) { # In stop-loss
pnls[i] <- 0 # Set PnLs = 0 if in stop-loss
minp <- min(minp, pricev[i]) # Update minimum price to current price
if (pricev[i] > ((1 + stopl)*minp)) { # Check for start-gain
insg <- TRUE # Is in start-gain?
insl <- FALSE # Is in stop-loss?
maxp <- pricev[i] # Reset trailing maximum price
} # end if
} else if (insg) { # In start-gain
maxp <- max(maxp, pricev[i]) # Update maximum price to current price
if (pricev[i] < ((1 - stopl)*maxp)) { # Check for stop-loss
insl <- TRUE # Is in stop-loss?
insg <- FALSE # Is in start-gain?
minp <- pricev[i] # Reset trailing minimum price
} # end if
} else { # Warmup period
# Update the maximum and minimum prices
maxp <- max(maxp, pricev[i])
minp <- min(minp, pricev[i])
# Update the stop-loss and start-gain indicators
insl <- (pricev[i] < ((1 - stopl)*maxp)) # Is in stop-loss?
insg <- (pricev[i] > ((1 + stopl)*minp)) # Is in start-gain?
} # end if
} # end for
return(pnls)
} # end sim_stopstart
# Simulate the stop-start strategy
pnls <- sim_stopstart(0.1)
# Combine the data
wealthv <- cbind(retp, pnls)
colnamev <- c("VTI", "Strategy")
colnames(wealthv) <- colnamev
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(wealthv, function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# dygraph plot the stop-loss strategy
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(cumsum(wealthv)[endd],
main="VTI Stop-Start Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Plot dygraph with shading
# Create colors for background shading
insl <- (pnls == 0) # Is in stop-loss?
indic <- (rutils::diffit(insl) != 0) # Indices of crosses
crossd <- c(datev[indic], datev[nrows]) # Dates of crosses
shadev <- ifelse(insl[indic] == 1, "antiquewhite", "lightgreen")
# Create dygraph object without plotting it
dyplot <- dygraphs::dygraph(cumsum(wealthv),
main="VTI Stop-Start Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Add shading to dygraph object
for (i in 1:NROW(shadev)) {
dyplot <- dyplot %>% dyShading(from=crossd[i], to=crossd[i+1], color=shadev[i])
} # end for
# Plot the dygraph object
dyplot
# Simulate multiple stop-loss strategies
stopv <- 0.01*(1:30)
pnlc <- sapply(stopv, function(stopl) {
sum(sim_stopstart(stopl))
}) # end sapply
stopl <- stopv[which.max(pnlc)]
# Plot cumulative pnls for stop-loss strategies
plot(x=stopv, y=pnlc,
main="Cumulative PnLs for Stop-Start Strategies",
xlab="stop-loss percent", ylab="cumulative pnl",
t="l", lwd=3, col="blue")
# Calculate the USO prices and returns
pricev <- na.omit(rutils::etfenv$prices$USO)
nrows <- NROW(pricev)
datev <- zoo::index(pricev)
retp <- rutils::diffit(log(pricev))
# Simulate multiple stop-start strategies
stopv <- 0.01*(1:30)
pnlc <- sapply(stopv, function(stopl) {
sum(sim_stopstart(stopl))
}) # end sapply
# Plot cumulative pnls for stop-start strategies
plot(x=stopv, y=pnlc,
main="Cumulative PnLs for USO Stop-Start Strategies",
xlab="stop-loss percent", ylab="cumulative pnl",
t="l", lwd=3, col="blue")
# Simulate optimal stop-start strategy for USO
stopl <- stopv[which.max(pnlc)]
pnls <- sim_stopstart(stopl)
# Combine the data
wealthv <- cbind(retp, pnls)
colnamev <- c("USO", "Strategy")
colnames(wealthv) <- colnamev
# Calculate the Sharpe and Sortino ratios
sqrt(252)*sapply(wealthv, function(x)
c(Sharpe=mean(x)/sd(x), Sortino=mean(x)/sd(x[x<0])))
# dygraph plot the stop-start strategy
endd <- rutils::calc_endpoints(wealthv, interval="weeks")
dygraphs::dygraph(cumsum(wealthv)[endd],
main="USO Stop-Start Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Plot dygraph with shading
# Create colors for background shading
insl <- (pnls == 0) # Is in stop-loss?
indic <- (rutils::diffit(insl) != 0) # Indices of crosses
crossd <- c(datev[indic], datev[nrows]) # Dates of crosses
shadev <- ifelse(insl[indic] == 1, "antiquewhite", "lightgreen")
# Create dygraph object without plotting it
dyplot <- dygraphs::dygraph(cumsum(wealthv),
main="USO Stop-Start Strategy") %>%
dyOptions(colors=c("blue", "red"), strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Add shading to dygraph object
for (i in 1:NROW(shadev)) {
dyplot <- dyplot %>% dyShading(from=crossd[i], to=crossd[i+1], color=shadev[i])
} # end for
# Plot the dygraph object
dyplot
# Extract the log VTI prices
ohlc <- log(rutils::etfenv$VTI)
closep <- quantmod::Cl(ohlc)
colnames(closep) <- "VTI"
nrows <- NROW(closep)
# Calculate the EMA weights
lookb <- 111
lambdaf <- 0.9
weightv <- lambdaf^(0:lookb)
weightv <- weightv/sum(weightv)
# Calculate the EMA prices as a convolution
pricema <- HighFreq::roll_sumw(closep, weightv=weightv)
pricev <- cbind(closep, pricema)
colnames(pricev) <- c("VTI", "VTI EMA")
# Dygraphs plot with custom line colors
colnamev <- colnames(pricev)
dygraphs::dygraph(pricev["2008/2009"], main="VTI EMA Prices") %>%
dySeries(name=colnamev[1], strokeWidth=1, col="blue") %>%
dySeries(name=colnamev[2], strokeWidth=2, col="red") %>%
dyLegend(show="always", width=300)
# Standard plot of EMA prices with custom line colors
x11(width=6, height=5)
plot_theme <- chart_theme()
colorv <- c("blue", "red")
plot_theme$col$line.col <- colorv
quantmod::chart_Series(pricev["2009"], theme=plot_theme,
lwd=2, name="VTI EMA Prices")
legend("topleft", legend=colnames(pricev), y.intersp=0.5,
inset=0.1, bg="white", lty=1, lwd=6, cex=0.8,
col=plot_theme$col$line.col, bty="n")
# Calculate the EMA prices recursively using C++ code
emar <- .Call(stats:::C_rfilter, closep, lambdaf, c(as.numeric(closep[1])/(1-lambdaf), double(NROW(closep))))[-1]
# Or R code
# emar <- filter(closep, filter=lambdaf, init=as.numeric(closep[1, 1])/(1-lambdaf), method="recursive")
emar <- (1-lambdaf)*emar
# Calculate the EMA prices recursively using RcppArmadillo
pricema <- HighFreq::run_mean(closep, lambda=lambdaf)
all.equal(drop(pricema), emar)
# Compare the speed of C++ code with RcppArmadillo
library(microbenchmark)
summary(microbenchmark(
Rcpp=HighFreq::run_mean(closep, lambda=lambdaf),
rfilter=.Call(stats:::C_rfilter, closep, lambdaf, c(as.numeric(closep[1])/(1-lambdaf), double(NROW(closep)))),
times=10))[, c(1, 4, 5)]
# Dygraphs plot with custom line colors
pricev <- cbind(closep, pricema)
colnames(pricev) <- c("VTI", "VTI EMA")
colnamev <- colnames(pricev)
dygraphs::dygraph(pricev["2008/2009"], main="Recursive VTI EMA Prices") %>%
dySeries(name=colnamev[1], strokeWidth=1, col="blue") %>%
dySeries(name=colnamev[2], strokeWidth=2, col="red") %>%
dyLegend(show="always", width=300)
# Standard plot of EMA prices with custom line colors
plot_theme <- chart_theme()
colorv <- c("blue", "red")
plot_theme$col$line.col <- colorv
quantmod::chart_Series(pricev["2009"], theme=plot_theme,
lwd=2, name="VTI EMA Prices")
legend("topleft", legend=colnames(pricev), y.intersp=0.5,
inset=0.1, bg="white", lty=1, lwd=6, cex=0.8,
col=plot_theme$col$line.col, bty="n")
# Calculate the EMA prices recursively using C++ code
lambdaf <- 0.984
pricema <- HighFreq::run_mean(closep, lambda=lambdaf)
pricev <- cbind(closep, pricema)
colnames(pricev) <- c("VTI", "VTI EMA")
colnamev <- colnames(pricev)
# Calculate the positions, either: -1, 0, or 1
indic <- sign(closep - pricema)
posv <- rutils::lagit(indic, lagg=1)
# Create colors for background shading
crossd <- (rutils::diffit(posv) != 0)
shadev <- posv[crossd]
crossd <- c(zoo::index(shadev), end(posv))
shadev <- ifelse(drop(zoo::coredata(shadev)) == 1, "lightgreen", "antiquewhite")
# Create dygraph object without plotting it
dyplot <- dygraphs::dygraph(pricev, main="VTI EMA Prices") %>%
dySeries(name=colnamev[1], strokeWidth=1, col="blue") %>%
dySeries(name=colnamev[2], strokeWidth=3, col="red") %>%
dyLegend(show="always", width=300)
# Add shading to dygraph object
for (i in 1:NROW(shadev)) {
dyplot <- dyplot %>% dyShading(from=crossd[i], to=crossd[i+1], color=shadev[i])
} # end for
# Plot the dygraph object
dyplot
# Standard plot of EMA prices with position shading
quantmod::chart_Series(pricev, theme=plot_theme,
lwd=2, name="VTI EMA Prices")
add_TA(posv > 0, on=-1, col="lightgreen", border="lightgreen")
add_TA(posv < 0, on=-1, col="lightgrey", border="lightgrey")
legend("topleft", legend=colnames(pricev),
inset=0.1, bg="white", lty=1, lwd=6, y.intersp=0.5,
col=plot_theme$col$line.col, bty="n")
# Calculate the daily profits and losses of crossover strategy
retp <- rutils::diffit(closep) # VTI returns
pnls <- retp*posv
colnames(pnls) <- "EMA"
wealthv <- cbind(retp, pnls)
colnames(wealthv) <- c("VTI", "EMA PnL")
# Annualized Sharpe ratio of crossover strategy
sqrt(252)*sapply(wealthv, function (x) mean(x)/sd(x))
# The crossover strategy has a negative correlation to VTI
cor(wealthv)[1, 2]
# Plot dygraph of crossover strategy wealth
# Create dygraph object without plotting it
colorv <- c("blue", "red")
dyplot <- dygraphs::dygraph(cumsum(wealthv), main="Performance of Crossover Strategy") %>%
dyOptions(colors=colorv, strokeWidth=2) %>%
dyLegend(show="always", width=300)
# Add shading to dygraph object
for (i in 1:NROW(shadev)) {
dyplot <- dyplot %>%
dyShading(from=crossd[i], to=crossd[i+1], color=shadev[i])
} # end for
# Plot the dygraph object
dyplot
# Standard plot of crossover strategy wealth
x11(width=6, height=5)
plot_theme <- chart_theme()
plot_theme$col$line.col <- colorv
quantmod::chart_Series(cumsum(wealthv), theme=plot_theme,
name="Performance of Crossover Strategy")
add_TA(posv > 0, on=-1, col="lightgreen", border="lightgreen")
add_TA(posv < 0, on=-1, col="lightgrey", border="lightgrey")
legend("top", legend=colnames(wealthv), y.intersp=0.5,
inset=0.05, bg="white", lty=1, lwd=6,
col=plot_theme$col$line.col, bty="n")
# Test EMA crossover market timing of VTI using Treynor-Mazuy test
desm <- cbind(pnls, retp, retp^2)
desm <- na.omit(desm)
colnames(desm) <- c("EMA", "VTI", "Treynor")
regmod <- lm(EMA ~ VTI + Treynor, data=desm)
summary(regmod)
# Plot residual scatterplot
resids <- (desm$EMA - regmod$coeff["VTI"]*retp)
resids <- regmod$residuals
plot.default(x=retp, y=resids, xlab="VTI", ylab="residuals")
title(main="Treynor-Mazuy Market Timing Test\n for EMA Crossover vs VTI", line=0.5)
# Plot fitted (predicted) response values
coefreg <- summary(regmod)$coeff
fitv <- regmod$fitted.values - coefreg["VTI", "Estimate"]*retp
tvalue <- round(coefreg["Treynor", "t value"], 2)
points.default(x=retp, y=fitv, pch=16, col="red")
text(x=0.0, y=0.8*max(resids), paste("Treynor test t-value =", tvalue))
# Determine trade dates right after EMA has crossed prices
indic <- sign(closep - pricema)
# Calculate the positions from lagged indicator
lagg <- 2
indic <- HighFreq::roll_sum(indic, lagg)
# Calculate the positions, either: -1, 0, or 1
posv <- rep(NA_integer_, nrows)
posv[1] <- 0
posv <- ifelse(indic == lagg, 1, posv)
posv <- ifelse(indic == (-lagg), -1, posv)
posv <- zoo::na.locf(posv, na.rm=FALSE)
posv <- xts::xts(posv, order.by=zoo::index(closep))
# Lag the positions to trade in next period
posv <- rutils::lagit(posv, lagg=1)
# Calculate the PnLs of lagged strategy
pnlslag <- retp*posv
colnames(pnlslag) <- "Lagged Strategy"
wealthv <- cbind(pnls, pnlslag)
colnames(wealthv) <- c("EMA", "Lagged")