options(digits=3, width=70)
library(IntroCompFinR)

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(PerformanceAnalytics)

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

library(tseries)
library(zoo)
library(boot)
library(corrplot)

## corrplot 0.84 loaded

setwd("~/Documents/UW")
getSymbols.warning4.0=FALSE
getSymbols.yahoo.warning=FALSE

# Gathering Data...
# 1. S&P 500 index (vfinx)
# 2. European stock index (veurx)
# 3. Emerging markets fund (veiex)
# 4. Long term bond index (vbltx)
# 5. Short term bond index (vbisx)
# 6. Pacific stock index (vpacx)
asset.names = c("vfinx","veurx","veiex","vbltx","vbisx","vpacx")
start.date = "2014-01-01"
end.date = "2019-01-31"
vfinx.prices = get.hist.quote(instrument="vfinx", start=start.date,
                             end=end.date, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")

## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.

## 
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
## 
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).

## time series ends   2019-01-01

veurx.prices = get.hist.quote(instrument="veurx", start=start.date,
                             end=end.date, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")

## time series ends   2019-01-01

veiex.prices = get.hist.quote(instrument="veiex", start=start.date,
                             end=end.date, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")

## time series ends   2019-01-01

vbltx.prices = get.hist.quote(instrument="vbltx", start=start.date,
                             end=end.date, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")

## time series ends   2019-01-01

vbisx.prices = get.hist.quote(instrument="vbisx", start=start.date,
                             end=end.date, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")

## time series ends   2019-01-01

vpacx.prices = get.hist.quote(instrument="vpacx", start=start.date,
                             end=end.date, quote="AdjClose",
                             provider="yahoo", origin="1970-01-01",
                             compression="m", retclass="zoo")

## time series ends   2019-01-01

# change time indices to class yearmon, which is most appropriate for monthly data
index(vfinx.prices) = as.yearmon(index(vfinx.prices))
index(veurx.prices) = as.yearmon(index(veurx.prices))
index(veiex.prices) = as.yearmon(index(veiex.prices))
index(vbltx.prices) = as.yearmon(index(vbltx.prices))
index(vbisx.prices) = as.yearmon(index(vbisx.prices))
index(vpacx.prices) = as.yearmon(index(vpacx.prices))
projectPrices.z = merge(vfinx.prices,veurx.prices,veiex.prices,vbltx.prices,
                        vbisx.prices,vpacx.prices)
colnames(projectPrices.z) = asset.names
# create data.frame for downloading
projectPrices.df = coredata(projectPrices.z)
rownames(projectPrices.df) = as.character(index(projectPrices.z))

RETURN CALCULATIONS & SAMPLE STATISTICS

# CC Returns
projectReturns.z = diff(log(projectPrices.z)) 
projectReturns.z

##             vfinx     veurx     veiex     vbltx     vbisx    vpacx
## Feb 2014  0.04460  0.069457  0.033624  0.014483  0.001805  0.03232
## Mar 2014  0.00395 -0.021567  0.035648  0.006347 -0.002948 -0.00265
## Apr 2014  0.01146  0.040062  0.007306  0.021022  0.003048 -0.00100
## May 2014  0.02308  0.007672  0.036492  0.024934  0.003704  0.03147
## Jun 2014  0.01611 -0.016025  0.020322  0.001246 -0.000950  0.02044
## Jul 2014 -0.00977 -0.026126  0.023724  0.003391 -0.001902  0.01388
## Aug 2014  0.03910  0.006136  0.035067  0.034532  0.002855 -0.00168
## Sep 2014 -0.01873 -0.045111 -0.084199 -0.025532 -0.001905 -0.05076
## Oct 2014  0.02839 -0.017453  0.032680  0.019997  0.004759  0.02755
## Nov 2014  0.02647  0.021118 -0.010043  0.018401  0.001900 -0.03423
## Dec 2014 -0.00792 -0.051169 -0.054940  0.015213 -0.004759 -0.03637
## Jan 2015 -0.02520  0.009055  0.010619  0.064648  0.011692  0.02905
## Feb 2015  0.05575  0.059002  0.034699 -0.036620 -0.004740  0.05483
## Mar 2015 -0.02086 -0.030095 -0.022223  0.004275  0.003704  0.00258
## Apr 2015  0.01428  0.046685  0.076230 -0.022899  0.000380  0.04097
## May 2015  0.01265  0.001645 -0.034373 -0.021783  0.001043 -0.00829
## Jun 2015 -0.02393 -0.048166 -0.032902 -0.035689 -0.001805 -0.02784
## Jul 2015  0.02498  0.039995 -0.061750  0.026375  0.001998 -0.00268
## Aug 2015 -0.06232 -0.072944 -0.097871 -0.015038 -0.000855 -0.07810
## Sep 2015 -0.03025 -0.045934 -0.046087  0.011688  0.004850 -0.04683
## Oct 2015  0.08582  0.063097  0.067306  0.005551 -0.001806  0.07577
## Nov 2015  0.00286 -0.012950 -0.033030 -0.004746 -0.001811 -0.00273
## Dec 2015 -0.02168 -0.028274 -0.029825 -0.012488 -0.003732 -0.01471
## Jan 2016 -0.04553 -0.056910 -0.057441  0.024578  0.009462 -0.04830
## Feb 2016 -0.00145 -0.028499 -0.008729  0.019654  0.002095 -0.03600
## Mar 2016  0.06048  0.059283  0.121272  0.029133  0.004841  0.07360
## Apr 2016  0.00883  0.033346  0.009060  0.013223  0.001138  0.01028
## May 2016  0.01767 -0.005613 -0.034590  0.000999 -0.001708  0.00468
## Jun 2016 -0.00222 -0.064313  0.043166  0.050127  0.009652 -0.00656
## Jul 2016  0.04085  0.056528  0.051636  0.025579  0.001132  0.06163
## Aug 2016  0.00129  0.007295  0.016607 -0.000974 -0.002548  0.00267
## Sep 2016 -0.00419  0.000765  0.000844 -0.010342  0.002174  0.02456
## Oct 2016 -0.01417 -0.029622  0.018104 -0.030545 -0.002650 -0.00726
## Nov 2016  0.03632 -0.023220 -0.044585 -0.064343 -0.009248 -0.01411
## Dec 2016  0.01345  0.042036 -0.005716  0.005385  0.000191 -0.01611
## Jan 2017  0.02463  0.033020  0.052285  0.008027  0.002399  0.05816
## Feb 2017  0.03886  0.006395  0.032248  0.017994  0.002204  0.02132
## Mar 2017 -0.00343  0.036086  0.020935 -0.006578  0.000191  0.01508
## Apr 2017  0.01451  0.045404  0.015184  0.014993  0.004115  0.00993
## May 2017  0.01384  0.046820  0.012107  0.018969  0.002197  0.02201
## Jun 2017  0.00183 -0.021441  0.001939  0.009563 -0.000574  0.01043
## Jul 2017  0.02448  0.043246  0.057480  0.003050  0.003153  0.03296
## Aug 2017  0.00293  0.001314  0.030437  0.020717  0.003244  0.00155
## Sep 2017  0.01551  0.027523 -0.019096 -0.010257 -0.002483  0.00773
## Oct 2017  0.02779  0.007741  0.035297  0.003756 -0.000574  0.04756
## Nov 2017  0.03013  0.001271  0.001774  0.005226 -0.003363  0.02044
## Dec 2017  0.00642  0.009485  0.030374  0.019034  0.000384  0.00217
## Jan 2018  0.06009  0.058606  0.084832 -0.020731 -0.005319  0.05872
## Feb 2018 -0.03769 -0.063088 -0.048425 -0.031611 -0.002431 -0.03647
## Mar 2018 -0.03014 -0.007973 -0.013410  0.012396  0.002337 -0.00645
## Apr 2018  0.00793  0.026221 -0.019730 -0.018746 -0.003316  0.00475
## May 2018  0.02366 -0.026990 -0.028673  0.009162  0.004494 -0.01008
## Jun 2018  0.00163 -0.029726 -0.052059 -0.002688  0.000584 -0.03385
## Jul 2018  0.04078  0.051713  0.037174  0.003130 -0.000391  0.01520
## Aug 2018  0.03195 -0.027328 -0.036140  0.006277  0.004592 -0.00895
## Sep 2018  0.00126 -0.002642 -0.024326 -0.013981 -0.002246  0.01487
## Oct 2018 -0.06669 -0.078744 -0.067643 -0.031457  0.000686 -0.09848
## Nov 2018  0.02006 -0.008293  0.043582  0.004221  0.002744  0.02102
## Dec 2018 -0.09984 -0.053926 -0.036530  0.037905  0.010531 -0.07386
## Jan 2019  0.08197  0.069690  0.088054  0.019090  0.004752  0.08463

# Simple Returns
projectReturnsSimple.z = exp(projectReturns.z) - 1
projectReturnsSimple.z

##             vfinx     veurx     veiex     vbltx     vbisx    vpacx
## Feb 2014  0.04561  0.071926  0.034195  0.014588  0.001807  0.03285
## Mar 2014  0.00395 -0.021336  0.036291  0.006367 -0.002943 -0.00265
## Apr 2014  0.01153  0.040875  0.007333  0.021245  0.003052 -0.00100
## May 2014  0.02335  0.007702  0.037166  0.025247  0.003711  0.03197
## Jun 2014  0.01624 -0.015897  0.020530  0.001247 -0.000949  0.02065
## Jul 2014 -0.00972 -0.025787  0.024007  0.003396 -0.001900  0.01398
## Aug 2014  0.03987  0.006155  0.035689  0.035136  0.002859 -0.00168
## Sep 2014 -0.01855 -0.044108 -0.080752 -0.025209 -0.001903 -0.04950
## Oct 2014  0.02879 -0.017302  0.033219  0.020198  0.004771  0.02794
## Nov 2014  0.02682  0.021343 -0.009992  0.018571  0.001902 -0.03365
## Dec 2014 -0.00789 -0.049882 -0.053458  0.015329 -0.004748 -0.03571
## Jan 2015 -0.02489  0.009096  0.010676  0.066784  0.011760  0.02947
## Feb 2015  0.05734  0.060777  0.035308 -0.035958 -0.004728  0.05636
## Mar 2015 -0.02064 -0.029647 -0.021978  0.004284  0.003711  0.00258
## Apr 2015  0.01439  0.047792  0.079211 -0.022639  0.000380  0.04182
## May 2015  0.01273  0.001647 -0.033789 -0.021547  0.001044 -0.00825
## Jun 2015 -0.02364 -0.047024 -0.032366 -0.035060 -0.001803 -0.02745
## Jul 2015  0.02529  0.040805 -0.059882  0.026726  0.002000 -0.00267
## Aug 2015 -0.06042 -0.070347 -0.093234 -0.014925 -0.000854 -0.07513
## Sep 2015 -0.02979 -0.044895 -0.045041  0.011757  0.004862 -0.04575
## Oct 2015  0.08961  0.065130  0.069622  0.005566 -0.001804  0.07872
## Nov 2015  0.00286 -0.012866 -0.032491 -0.004735 -0.001810 -0.00273
## Dec 2015 -0.02144 -0.027878 -0.029384 -0.012410 -0.003725 -0.01460
## Jan 2016 -0.04451 -0.055321 -0.055823  0.024883  0.009507 -0.04715
## Feb 2016 -0.00145 -0.028097 -0.008691  0.019848  0.002097 -0.03536
## Mar 2016  0.06235  0.061075  0.128932  0.029562  0.004853  0.07637
## Apr 2016  0.00887  0.033909  0.009101  0.013311  0.001139  0.01034
## May 2016  0.01783 -0.005597 -0.033998  0.000999 -0.001706  0.00469
## Jun 2016 -0.00222 -0.062289  0.044111  0.051405  0.009699 -0.00654
## Jul 2016  0.04169  0.058156  0.052993  0.025909  0.001133  0.06357
## Aug 2016  0.00130  0.007322  0.016745 -0.000974 -0.002545  0.00267
## Sep 2016 -0.00418  0.000765  0.000845 -0.010289  0.002177  0.02487
## Oct 2016 -0.01407 -0.029188  0.018269 -0.030083 -0.002647 -0.00724
## Nov 2016  0.03699 -0.022952 -0.043606 -0.062316 -0.009206 -0.01401
## Dec 2016  0.01354  0.042933 -0.005699  0.005400  0.000191 -0.01599
## Jan 2017  0.02493  0.033571  0.053676  0.008060  0.002402  0.05989
## Feb 2017  0.03963  0.006415  0.032773  0.018157  0.002207  0.02155
## Mar 2017 -0.00343  0.036745  0.021155 -0.006556  0.000191  0.01519
## Apr 2017  0.01461  0.046451  0.015300  0.015106  0.004124  0.00998
## May 2017  0.01394  0.047933  0.012181  0.019150  0.002199  0.02226
## Jun 2017  0.00184 -0.021213  0.001941  0.009608 -0.000573  0.01048
## Jul 2017  0.02478  0.044195  0.059164  0.003055  0.003158  0.03351
## Aug 2017  0.00293  0.001315  0.030905  0.020933  0.003249  0.00156
## Sep 2017  0.01563  0.027905 -0.018915 -0.010204 -0.002480  0.00776
## Oct 2017  0.02818  0.007771  0.035927  0.003763 -0.000574  0.04871
## Nov 2017  0.03059  0.001272  0.001776  0.005240 -0.003357  0.02065
## Dec 2017  0.00644  0.009530  0.030840  0.019217  0.000384  0.00217
## Jan 2018  0.06193  0.060358  0.088535 -0.020518 -0.005305  0.06048
## Feb 2018 -0.03698 -0.061139 -0.047272 -0.031117 -0.002428 -0.03581
## Mar 2018 -0.02969 -0.007942 -0.013320  0.012473  0.002340 -0.00643
## Apr 2018  0.00796  0.026568 -0.019537 -0.018572 -0.003311  0.00476
## May 2018  0.02395 -0.026629 -0.028266  0.009204  0.004504 -0.01003
## Jun 2018  0.00164 -0.029289 -0.050727 -0.002685  0.000584 -0.03329
## Jul 2018  0.04162  0.053074  0.037873  0.003135 -0.000391  0.01532
## Aug 2018  0.03246 -0.026958 -0.035495  0.006297  0.004602 -0.00891
## Sep 2018  0.00127 -0.002639 -0.024033 -0.013884 -0.002244  0.01498
## Oct 2018 -0.06452 -0.075724 -0.065406 -0.030968  0.000686 -0.09379
## Nov 2018  0.02027 -0.008258  0.044546  0.004230  0.002747  0.02124
## Dec 2018 -0.09502 -0.052498 -0.035871  0.038633  0.010586 -0.07120
## Jan 2019  0.08543  0.072176  0.092047  0.019273  0.004763  0.08831

# Sorting into a data frame...
projectReturns.df = as.data.frame(coredata(projectReturns.z))
rownames(projectReturns.df) = as.character(index(projectReturns.z))
projectReturnsSimple.df = as.data.frame(coredata(projectReturnsSimple.z))
rownames(projectReturnsSimple.df) = as.character(index(projectReturnsSimple.z))

Time Plots of Monthly Prices and Continuously Compounded Returns

# Plotting Data...
my.panel <- function(...) {
  lines(...)
  abline(h=0)
}
plot(projectPrices.z, col="blue", lwd=2)

plot(projectReturns.z, panel=my.panel, col="blue", lwd=2)

VFINX, VBLTX, and VBISX experience similar, positive growths in prices and experience less volatility overall. VEURX, VEIEX, and VPACX experience a positive growth from 2014 to 2019 as well but it appears more volatile, experience primarily negative growth from 2014 to 2016 until price growth turns positive. A likely explanation can be the 2015-2016 stock market selloff, when a period of decline in the value of stock prices occurred globally. This time period included the Chinese stock market turbulence where the SSE Composite Index fell 43% over two months, triggering the devaluation of the yuan, a slowing of growth in China’s GDP, and a large sell-off of assets from investors globally.

All six of the assets experience large and sharp increases or decreases in return rates at the same time. The smaller and less sharp increases or decreases in return rates are not always in the same negative or positive direction.

Assets VPACX, VEIEX, VFINX, and VEURX are the more volatile assets and experience the largest increases and decreases in return rates (about +/-10%). VBLTX is relatively less volatile and experiences return rates between +/- 6%. VBISX experiences the lowest amount of volatility as their return rates fluctuate between +/-1%.

Overall they follow very similar trends, all experiencing a sharp increases in return rates in beginning-2015 and 2019, as well as sharp decreases in mid-2015 and end-2018.

Equity Curves

# $1 Growth over 5 Years...
chart.CumReturns(projectReturnsSimple.z, wealth.index=TRUE, legend.loc="topleft", 
                 lwd=2, main="growth of $1") 

Upon examining the equity curves, VFINX gives the highest future value of about $250. This is to be expected since VFINX tracks returns on investments for the S&P 500 index which has experienced a substantial growth since 2014.

Four-panel Diagnostic Plots

ret.mat = coredata(projectReturns.z)
# Four-panel Plot for Vangaurd Assets
fourPanelPlot(projectReturns.z[, "vfinx", drop=FALSE])

fourPanelPlot(projectReturns.z[, "veurx", drop=FALSE])

fourPanelPlot(projectReturns.z[, "veiex", drop=FALSE])

fourPanelPlot(projectReturns.z[, "vbltx", drop=FALSE])

fourPanelPlot(projectReturns.z[, "vbisx", drop=FALSE])

fourPanelPlot(projectReturns.z[, "vpacx", drop=FALSE])

# Just Boxplots...
boxplot(ret.mat, main="Vanguard Returns", col="blue")

None of the returns appear to be perfectly normally distributed. The Q-Q plots show that all assets have fat tails. VEURX appears to fit a normal distribution the least while VBlTX and VEIEX fit one the most. The distribution plots show that VEIEX and VBISX experience a positive skew. The box-plots show that VFINX, VBLTX, VBISX, and VPACX all have outliers. Although the ACF plots show that all assets have a weak time-dependence, VBLTC and VEURX show more time-sensitivity in comparison to the other assets.

Univariate Descritive Statistics

# Descriptive Statistics (CER MODEL)
muhat.vals = colMeans(projectReturns.z)
var.vals = apply(projectReturns.z, 2, var)
sd.vals = apply(projectReturns.z, 2, sd)
skew.vals = apply(projectReturns.z, 2, skewness)
ekurt.vals = apply(projectReturns.z, 2, kurtosis)
# Empirical Quantiles
q.vals = apply(projectReturns.z, 2, quantile, prob=c(0.01,0.05))
# Organizing Results into a Table
stats.mat = rbind(muhat.vals,
                  var.vals,
                  sd.vals,
                  skew.vals,
                  ekurt.vals,
                  q.vals)
rownames(stats.mat) = c("Mean","Variance", "Std Dev", "Skewness", 
                        "Excess Kurtosis", "1% Quantile", 
                        "5% Quantile")
stats.mat

##                    vfinx    veurx    veiex     vbltx     vbisx
## Mean             0.00854  0.00139  0.00356  0.003954  9.49e-04
## Variance         0.00112  0.00157  0.00201  0.000505  1.53e-05
## Std Dev          0.03342  0.03962  0.04487  0.022477  3.91e-03
## Skewness        -0.50676 -0.04748  0.08554 -0.311537  4.29e-01
## Excess Kurtosis  1.27721 -0.95176 -0.25522  0.811702  7.76e-01
## 1% Quantile     -0.08028 -0.07532 -0.08980 -0.047986 -6.93e-03
## 5% Quantile     -0.04637 -0.06315 -0.06204 -0.031815 -4.74e-03
##                    vpacx
## Mean             0.00417
## Variance         0.00136
## Std Dev          0.03687
## Skewness        -0.27553
## Excess Kurtosis  0.41607
## 1% Quantile     -0.08646
## 5% Quantile     -0.05192

VFINX has the highest average return while VBISX has the lowest. VEIEX has the highest standard deviation while VBISX has the lowest. VBlTX and VEIEX look the most normally distributed while VEURX looks the least.

Sharpe’s Slope [table in excel]

# Set risk-free rate of 0.0004167 per month or 0.5% CC annual rate...
rf = 0.005/12
SharpeRatios = (muhat.vals - rf)/sd.vals
SharpeRatios

##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## 0.2432 0.0246 0.0701 0.1574 0.1362 0.1019

# Bootstrap Estimated Standard Errors...
# sharpe ratio bootstrap function...
sharpeRatio.boot = function(x, idx, risk.free) {
  muhat = mean(x[idx])
  sigmahat = sd(x[idx])
  sharpeRatio = (muhat - risk.free)/sigmahat
  sharpeRatio
}
# VFINX bootstrap
sharpe.vfinx.boot = boot(ret.mat[, "vfinx"], 
                         statistic=sharpeRatio.boot, R=999, risk.free=rf)
sharpe.vfinx.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vfinx"], statistic = sharpeRatio.boot, 
##     R = 999, risk.free = rf)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    0.243 0.00836       0.149

# VEURX bootstrap
sharpe.veurx.boot = boot(ret.mat[, "veurx"], 
                         statistic=sharpeRatio.boot, R=999, risk.free=rf)
sharpe.veurx.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "veurx"], statistic = sharpeRatio.boot, 
##     R = 999, risk.free = rf)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*   0.0246 0.00289       0.129

# VEIEX bootstrap
sharpe.veiex.boot = boot(ret.mat[, "veiex"], 
                         statistic=sharpeRatio.boot, R=999, risk.free=rf)
sharpe.veiex.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "veiex"], statistic = sharpeRatio.boot, 
##     R = 999, risk.free = rf)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*   0.0701 0.00784       0.134

# VBLTX bootstrap
sharpe.vbltx.boot = boot(ret.mat[, "vbltx"], 
                         statistic=sharpeRatio.boot, R=999, risk.free=rf)
sharpe.vbltx.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vbltx"], statistic = sharpeRatio.boot, 
##     R = 999, risk.free = rf)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    0.157 0.00788       0.134

# VBISX bootstrap
sharpe.vbisx.boot = boot(ret.mat[, "vbisx"], 
                         statistic=sharpeRatio.boot, R=999, risk.free=rf)
sharpe.vbisx.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vbisx"], statistic = sharpeRatio.boot, 
##     R = 999, risk.free = rf)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    0.136 -0.0058       0.127

# VPACX bootstrap
sharpe.vpacx.boot = boot(ret.mat[, "vpacx"], 
                         statistic=sharpeRatio.boot, R=999, risk.free=rf)
sharpe.vpacx.boot

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vpacx"], statistic = sharpeRatio.boot, 
##     R = 999, risk.free = rf)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    0.102 0.00926       0.133

VFINX has the highest Sharpe Ratio. Given the small standard errors, the Sharpe ratios appear to be well-estimated.

Standard Errors and 95% Confidence Intervals for Mean and SD [table in excel]

cov.mat = var(projectReturns.z)
cor.mat = cov2cor(cov.mat)
covhat.vals = cov.mat[lower.tri(cov.mat)]
rhohat.vals = cor.mat[lower.tri(cor.mat)]
names(covhat.vals) <- names(rhohat.vals) <- 
  c("vfinx,veurx","vfinx,veiex","vfinx,vbltx", "vfinx,vbisx", "vfinx,vpacx",
    "veurx,veiex", "veurx,vbltx", "veurx,vbisx", "veurx,vpacx",
    "veiex,vbltx", "veiex,vbisx", "veiex,vpacx",
    "vbltx,vbisx", "vbltx,vpacx",
    "vbisx,vpacx")
# Estimates for CER Model Parameters µi, σ2i, σi, σij, and ρij
muhat.vals = apply(ret.mat, 2, mean)
muhat.vals

##    vfinx    veurx    veiex    vbltx    vbisx    vpacx 
## 0.008545 0.001393 0.003563 0.003954 0.000949 0.004175

sigma2hat.vals = apply(ret.mat, 2, var)
sigma2hat.vals

##    vfinx    veurx    veiex    vbltx    vbisx    vpacx 
## 1.12e-03 1.57e-03 2.01e-03 5.05e-04 1.53e-05 1.36e-03

sigmahat.vals = apply(ret.mat, 2, sd)
sigmahat.vals

##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
## 0.03342 0.03962 0.04487 0.02248 0.00391 0.03687

cor.mat = cor(ret.mat)
cor.mat

##         vfinx   veurx veiex  vbltx   vbisx   vpacx
## vfinx  1.0000  0.7558 0.662 0.0509 -0.2129  0.8039
## veurx  0.7558  1.0000 0.680 0.1483 -0.0588  0.7875
## veiex  0.6615  0.6800 1.000 0.2902  0.1195  0.8414
## vbltx  0.0509  0.1483 0.290 1.0000  0.7679  0.1629
## vbisx -0.2129 -0.0588 0.119 0.7679  1.0000 -0.0423
## vpacx  0.8039  0.7875 0.841 0.1629 -0.0423  1.0000

rhohat.vals = cor.mat[lower.tri(cor.mat)]
names(rhohat.vals) <- 
c("VFINX,VEURX","VFINX,VEIEX","VFINX,VBLTX","VFINX,VBISX","VFINX,VPACX","VEURX,VEIEX","VEURX,VBLTX","VEURX,VBISX","VEURX,VPACX","VEIEX,VBLTX","VEIEX,VBISX","VEIEX,VPACX","VBLTX,VBISX","VBLTX,VPACX","VBISX,VPACX")
rhohat.vals

## VFINX,VEURX VFINX,VEIEX VFINX,VBLTX VFINX,VBISX VFINX,VPACX 
##      0.7558      0.6615      0.0509     -0.2129      0.8039 
## VEURX,VEIEX VEURX,VBLTX VEURX,VBISX VEURX,VPACX VEIEX,VBLTX 
##      0.6800      0.1483     -0.0588      0.7875      0.2902 
## VEIEX,VBISX VEIEX,VPACX VBLTX,VBISX VBLTX,VPACX VBISX,VPACX 
##      0.1195      0.8414      0.7679      0.1629     -0.0423

cbind(muhat.vals,sigma2hat.vals,sigmahat.vals)

##       muhat.vals sigma2hat.vals sigmahat.vals
## vfinx   0.008545       1.12e-03       0.03342
## veurx   0.001393       1.57e-03       0.03962
## veiex   0.003563       2.01e-03       0.04487
## vbltx   0.003954       5.05e-04       0.02248
## vbisx   0.000949       1.53e-05       0.00391
## vpacx   0.004175       1.36e-03       0.03687

cbind(covhat.vals,rhohat.vals)

##             covhat.vals rhohat.vals
## vfinx,veurx    1.00e-03      0.7558
## vfinx,veiex    9.92e-04      0.6615
## vfinx,vbltx    3.82e-05      0.0509
## vfinx,vbisx   -2.78e-05     -0.2129
## vfinx,vpacx    9.90e-04      0.8039
## veurx,veiex    1.21e-03      0.6800
## veurx,vbltx    1.32e-04      0.1483
## veurx,vbisx   -9.11e-06     -0.0588
## veurx,vpacx    1.15e-03      0.7875
## veiex,vbltx    2.93e-04      0.2902
## veiex,vbisx    2.10e-05      0.1195
## veiex,vpacx    1.39e-03      0.8414
## vbltx,vbisx    6.75e-05      0.7679
## vbltx,vpacx    1.35e-04      0.1629
## vbisx,vpacx   -6.10e-06     -0.0423

# Estimated Standard Errors for µi, σ2i, σi, and ρij
# SE(muhat)
nobs = nrow(ret.mat)
se.muhat = sigmahat.vals/sqrt(nobs)
rbind(muhat.vals,se.muhat)

##              vfinx   veurx   veiex   vbltx    vbisx   vpacx
## muhat.vals 0.00854 0.00139 0.00356 0.00395 0.000949 0.00417
## se.muhat   0.00431 0.00511 0.00579 0.00290 0.000505 0.00476

# SE(sigmahat) and SE(sigma2hat)
se.sigma2hat = sigma2hat.vals/sqrt(nobs/2)
se.sigmahat = sigmahat.vals/sqrt(2*nobs)
rbind(sigma2hat.vals,se.sigma2hat)

##                   vfinx    veurx    veiex    vbltx    vbisx    vpacx
## sigma2hat.vals 0.001117 0.001570 0.002013 5.05e-04 1.53e-05 0.001359
## se.sigma2hat   0.000204 0.000287 0.000368 9.22e-05 2.79e-06 0.000248

rbind(sigmahat.vals,se.sigmahat)

##                 vfinx   veurx  veiex   vbltx    vbisx   vpacx
## sigmahat.vals 0.03342 0.03962 0.0449 0.02248 0.003910 0.03687
## se.sigmahat   0.00305 0.00362 0.0041 0.00205 0.000357 0.00337

# SE(rhohat)
se.rhohat = (1-rhohat.vals^2)/sqrt(nobs)
rbind(rhohat.vals,se.rhohat)

##             VFINX,VEURX VFINX,VEIEX VFINX,VBLTX VFINX,VBISX
## rhohat.vals      0.7558      0.6615      0.0509      -0.213
## se.rhohat        0.0554      0.0726      0.1288       0.123
##             VFINX,VPACX VEURX,VEIEX VEURX,VBLTX VEURX,VBISX
## rhohat.vals      0.8039      0.6800       0.148     -0.0588
## se.rhohat        0.0457      0.0694       0.126      0.1287
##             VEURX,VPACX VEIEX,VBLTX VEIEX,VBISX VEIEX,VPACX
## rhohat.vals       0.787       0.290       0.119      0.8414
## se.rhohat         0.049       0.118       0.127      0.0377
##             VBLTX,VBISX VBLTX,VPACX VBISX,VPACX
## rhohat.vals       0.768       0.163     -0.0423
## se.rhohat         0.053       0.126      0.1289

# SE(muhat)
rbind(muhat.vals,se.muhat)

##              vfinx   veurx   veiex   vbltx    vbisx   vpacx
## muhat.vals 0.00854 0.00139 0.00356 0.00395 0.000949 0.00417
## se.muhat   0.00431 0.00511 0.00579 0.00290 0.000505 0.00476

# SE(sigmahat)
rbind(sigmahat.vals,se.sigmahat)

##                 vfinx   veurx  veiex   vbltx    vbisx   vpacx
## sigmahat.vals 0.03342 0.03962 0.0449 0.02248 0.003910 0.03687
## se.sigmahat   0.00305 0.00362 0.0041 0.00205 0.000357 0.00337

# 95% Confidence Intervals
# For mean...
mu.lower.95 = muhat.vals - 2*se.muhat
mu.upper.95 = muhat.vals + 2*se.muhat
cbind(mu.lower.95,mu.upper.95)

##       mu.lower.95 mu.upper.95
## vfinx   -8.34e-05     0.01717
## veurx   -8.84e-03     0.01162
## veiex   -8.02e-03     0.01515
## vbltx   -1.85e-03     0.00976
## vbisx   -6.02e-05     0.00196
## vpacx   -5.34e-03     0.01369

# For SD...
sigma.lower.95 = sigmahat.vals - 2*se.sigmahat
sigma.upper.95 = sigmahat.vals + 2*se.sigmahat
cbind(sigma.lower.95,sigma.upper.95)

##       sigma.lower.95 sigma.upper.95
## vfinx         0.0273        0.03952
## veurx         0.0324        0.04685
## veiex         0.0367        0.05306
## vbltx         0.0184        0.02658
## vbisx         0.0032        0.00462
## vpacx         0.0301        0.04360

The 95% confidence intervals for the estimated mean values are negative and positive values, so there is uncertainty about whether the true values of the mean returns are positive or negative.The 95% confidence intervals for the estimated standard deviation values and mean values are significantly narrow for all six assets, however. Overall, I would say that the estimations for the means and standard deviations are precise. However, the standard deviations are estimated more precisely than the means.

Sharpe’s Slope Annualized [table in excel]

# Annualize Monthly Statistics...
muhat.annual = 12*muhat.vals
sd.annual = sqrt(12)*sd.vals # <--- this is sqrt(12*var)
muhat.annual

##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## 0.1025 0.0167 0.0428 0.0474 0.0114 0.0501

sd.annual

##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## 0.1158 0.1372 0.1554 0.0779 0.0135 0.1277

# Annualized Sharpe Ratios
SharpeRatios.annual = (muhat.annual - rf)/sd.annual
SharpeRatios.annual

## vfinx veurx veiex vbltx vbisx vpacx 
## 0.882 0.119 0.272 0.604 0.810 0.389

# CC FV $1 over 5 Years in excel

The annualized mean returns are smaller for each asset in comparison to the monthly mean returns. The annualized standard deviations are larger than the monthly standard deviation values for each asset. VFINX still has the highest Sharpe slope and remains the preferred asset.

Pair-wise Scatterplots & Sample Covariance Matrix [table in excel]

# Matrix of Return Data and Pairwise Scatterplots
pairs(ret.mat, col="blue")

# how to write covariance matrix, expected returns, sd values and quantiles to files for importing into Excel (already ran)
# write.csv(projectPrices.df, file=paste("projectPrices.csv", sep=""))
# write.csv(projectReturns.df, file=paste("projectReturns.csv", sep=""))
# write.csv(projectReturnsSimple.df, file=paste("projectReturns.csv", sep=""))
# write.csv(cov.mat, file=paste("covmat.csv", sep=""))
# write.csv(cor.mat, file=paste("cormat.csv", sep=""))
# write.csv(muhat.vals, file=paste("muhat.csv", sep=""))
# write.csv(sd.vals, file=paste("sd.csv", sep=""))
# write.csv(t(q.vals), file=paste("q.csv", sep=""))

VFINX and VEURX , VFINX and VPACX, and VEIEX and VPACX show strong positive relationships. All other assets appear to have very weak correlated relationships. The sample covariance matrix shows negative relationships between VBISX and VFINX, VBISX and VEURX, VBISX and VPACX. All other assets have positive relationships, and many have very weak relationships.

Sample Correlation Matrix [table in excel]

# Plot
corrplot(cor.mat, method="ellipse")

VPACX and VEURX, VPACX and VEIEX, and VBLTX and VBISX are the most highly correlated assets. VFINX and VBLTX are the least correlated assets. Since the correlations vary incredibly with each asset, risk should be reduced with diversification of these assets.

VALUE-AT-RISK CALCULATIONS

1% and 5% value-at-risk of the $100,000 investment over one month

# Normal VaR Function
Value.at.Risk = function(x, p=0.05, w=100000, method=c("normal", "empirical")) {
    method=method[1]
  x = as.matrix(x)
  if (method == "normal") {
      q = apply(x, 2, mean) + apply(x, 2, sd)*qnorm(p)
  } else {    
    q = apply(x, 2, quantile, p)
  }
    VaR = (exp(q) - 1)*w
    VaR
}
# 5% and 1% Normal VaR for All Assets
VaR.normal.05 = Value.at.Risk(ret.mat, p=0.05, method="normal")
VaR.normal.05

## vfinx veurx veiex vbltx vbisx vpacx 
## -4536 -6178 -6783 -3248  -547 -5490

VaR.normal.01 = Value.at.Risk(ret.mat, p=0.01, method="normal")
VaR.normal.01

## vfinx veurx veiex vbltx vbisx vpacx 
## -6685 -8678 -9590 -4719  -811 -7835

Bootstrap Estimated Standard Errors & 95% Confidence Intervals for 5% and 1% VaR Estimates

VFINX

# Function for 5%...
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.05.boot.vfinx = boot(ret.mat[, "vfinx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.05.boot.vfinx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vfinx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -4536     115         797

boot.ci(VaR.05.boot.vfinx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.05.boot.vfinx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-6214, -3088 )   (-5975, -2965 )  
## Calculations and Intervals on Original Scale

# Function for 1%
ValueAtRisk.boot = function(x, idx, p=0.01, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.01.boot.vfinx = boot(ret.mat[, "vfinx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.01.boot.vfinx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vfinx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -6685     125         995

boot.ci(VaR.01.boot.vfinx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.01.boot.vfinx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-8760, -4860 )   (-8567, -4672 )  
## Calculations and Intervals on Original Scale

VEURX

# Function for 5%...
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.05.boot.veurx = boot(ret.mat[, "veurx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.05.boot.veurx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "veurx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -6178    75.6         647

boot.ci(VaR.05.boot.veurx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.05.boot.veurx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-7521, -4987 )   (-7313, -4722 )  
## Calculations and Intervals on Original Scale

# Function for 1%
ValueAtRisk.boot = function(x, idx, p=0.01, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.01.boot.veurx = boot(ret.mat[, "veurx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.01.boot.veurx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "veurx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -8678     139         711

boot.ci(VaR.01.boot.veurx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.01.boot.veurx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-10212,  -7423 )   ( -9909,  -7135 )  
## Calculations and Intervals on Original Scale

VEIEX

# Function for 5%...
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.05.boot.veiex = boot(ret.mat[, "veiex"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.05.boot.veiex

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "veiex"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -6783     124         786

boot.ci(VaR.05.boot.veiex, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.05.boot.veiex, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-8447, -5367 )   (-8175, -5010 )  
## Calculations and Intervals on Original Scale

# Function for 1%
ValueAtRisk.boot = function(x, idx, p=0.01, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.01.boot.veiex = boot(ret.mat[, "veiex"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.01.boot.veiex

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "veiex"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -9590    86.2         950

boot.ci(VaR.01.boot.veiex, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.01.boot.veiex, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-11538,  -7815 )   (-11321,  -7527 )  
## Calculations and Intervals on Original Scale

VBLTX

# Function for 5%...
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.05.boot.vbltx = boot(ret.mat[, "vbltx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.05.boot.vbltx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vbltx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -3248    65.2         509

boot.ci(VaR.05.boot.vbltx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.05.boot.vbltx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-4312, -2314 )   (-4221, -2207 )  
## Calculations and Intervals on Original Scale

# Function for 1%
ValueAtRisk.boot = function(x, idx, p=0.01, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.01.boot.vbltx = boot(ret.mat[, "vbltx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.01.boot.vbltx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vbltx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -4719    59.7         672

boot.ci(VaR.01.boot.vbltx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.01.boot.vbltx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-6095, -3462 )   (-5979, -3374 )  
## Calculations and Intervals on Original Scale

VBISX

# Function for 5%...
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.05.boot.vbisx = boot(ret.mat[, "vbisx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.05.boot.vbisx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vbisx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     -547    7.78        74.1

boot.ci(VaR.05.boot.vbisx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.05.boot.vbisx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-700, -409 )   (-686, -398 )  
## Calculations and Intervals on Original Scale

# Function for 1%
ValueAtRisk.boot = function(x, idx, p=0.01, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.01.boot.vbisx = boot(ret.mat[, "vbisx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.01.boot.vbisx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vbisx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     -811    15.9        97.2

boot.ci(VaR.01.boot.vbisx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.01.boot.vbisx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-1018,  -637 )   ( -989,  -615 )  
## Calculations and Intervals on Original Scale

VPACX

# Function for 5%...
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.05.boot.vpacx = boot(ret.mat[, "vpacx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.05.boot.vpacx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vpacx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -5490    59.1         784

boot.ci(VaR.05.boot.vpacx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.05.boot.vpacx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-7086, -4013 )   (-6975, -3791 )  
## Calculations and Intervals on Original Scale

# Function for 1%
ValueAtRisk.boot = function(x, idx, p=0.01, w=100000) {

    q = mean(x[idx]) + sd(x[idx])*qnorm(p)
    VaR = (exp(q) - 1)*w
    VaR
}
VaR.01.boot.vpacx = boot(ret.mat[, "vpacx"], 
                         statistic=ValueAtRisk.boot, R=999)
VaR.01.boot.vpacx

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = ret.mat[, "vpacx"], statistic = ValueAtRisk.boot, 
##     R = 999)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    -7835    92.6         947

boot.ci(VaR.01.boot.vpacx, conf = 0.95, type = c("norm","perc"))

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 999 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = VaR.01.boot.vpacx, conf = 0.95, type = c("norm", 
##     "perc"))
## 
## Intervals : 
## Level      Normal             Percentile     
## 95%   (-9784, -6071 )   (-9640, -5905 )  
## Calculations and Intervals on Original Scale

VEIEX has the highest 5% and 1% Value-at-Risk values. VBISX has the lowest 5% and 1% Value-at-Risk values. VFINX has the largest standard error for both the 5% and 1% Value-at-Risk values. VBISX has the smallest standard error for both the 5% and 1% Value-at-Risk values. The confidence intervals for the 5% VaR estimates are considerably wide, which is cause for concern about the accuracy of the estimations. However, the confidence intervals for the 1% VaR estimates are less wide, therefor indicating that these estimations are more accurate.

1% and 5% value-at-risk of the $100,000 investment over one year [table in excel]

Value.at.Risk = function(x, p=0.05, w=100000, method=c("normal", "empirical")) {
    method=method[1]
  x = as.matrix(x)
  if (method == "normal") {
      q = apply(x, 2, mean)*12 + apply(x, 2, sd)*sqrt(12)*qnorm(p)
  } else {    
    q = apply(x, 2, quantile, p)
  }
    VaR = (exp(q) - 1)*w
    VaR
}
# 5% and 1% Normal VaR for All Assets
VaR.normal.05 = Value.at.Risk(ret.mat, p=0.05, method="normal")
VaR.normal.05

##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
##  -8412 -18863 -19177  -7746  -1083 -14783

VaR.normal.01 = Value.at.Risk(ret.mat, p=0.01, method="normal")
VaR.normal.01

##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## -15359 -26108 -27300 -12513  -1992 -21887

Empirical/Historical VaR

# empirical VaR
VaR.empirical.05 = Value.at.Risk(ret.mat, p=0.05, method="empirical")
VaR.empirical.05

## vfinx veurx veiex vbltx vbisx vpacx 
## -4531 -6120 -6016 -3131  -473 -5059

VaR.empirical.01 = Value.at.Risk(ret.mat, p=0.01, method="empirical")
VaR.empirical.01

## vfinx veurx veiex vbltx vbisx vpacx 
## -7714 -7256 -8589 -4685  -691 -8283

The 5% VaR values using an empirical distribution are slightly smaller than those values calculated with a normal distribution. The 1% VaR values using an an empirical distribution vary from those values calculated with a normal distribution.

PORTFOLIO THEORY

Global Minimum Variance Portfolio [table in excel]

gmin.port <- globalMin.portfolio(muhat.vals, cov.mat)
summary(gmin.port, risk.free=rf)

## Call:
## globalMin.portfolio(er = muhat.vals, cov.mat = cov.mat)
## 
## Portfolio expected return:     0.000897 
## Portfolio standard deviation:  0.00251 
## Portfolio Sharpe Ratio:        0.191 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.0484  0.0003 -0.0151 -0.1291  1.0932  0.0024

plot(gmin.port)

In the global minimum variance portfolio, VEIEX and VBLTX have negative weights and thus require shorting the assets.

Annual Sharpe Ratio [table in excel]

SharpeRatio.glblmin.annual = ((0.000897*12) - rf)/(0.00251*sqrt(12))
SharpeRatio.glblmin.annual

## [1] 1.19

The yearly expected return and standard deviation for the global minimum variance portfolio are less than that of each asset individually. The annual Sharpe ratio, however, is larger than each of the assets’ individual annual Sharpe ratios. This gives the indication that this portfolio is preferrable to each of the assets individually.

One-Month VaR [table in excel]

The Value-at-Risk value for the global minimum variance portfolio is significantly smaller than the VaR values for each individual asset.

Global Minimum Variance Portfolio w/ No Shorts [table in excel]

gmin.port.ns <- globalMin.portfolio(muhat.vals, cov.mat, shorts=FALSE)
summary(gmin.port.ns, risk.free=rf)

## Call:
## globalMin.portfolio(er = muhat.vals, cov.mat = cov.mat, shorts = FALSE)
## 
## Portfolio expected return:     0.00123 
## Portfolio standard deviation:  0.0037 
## Portfolio Sharpe Ratio:        0.218 
## Portfolio weights:
##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## 0.0363 0.0000 0.0000 0.0000 0.9637 0.0000

plot(gmin.port.ns)

Annual Sharpe Ratio [table in excel]

The global minimum variance portfolio without short sales has a higher mean, standard deviation, and sharpe ratio than those of the global minimum variance portfolio with short sales.

One-Month VaR [table in excel]

The 1% and 5% VaR values are larger in the no-shorts global minimum variance portfolio than those of those of the global minimum variance portfolio with short sales.

Efficient Portfolio Frontier w/ Short Sales

ef <- efficient.frontier(muhat.vals, cov.mat, alpha.min=-1, 
                         alpha.max=1.5, nport=20)
# Efficient portfolio with target return equal to max returns...
target.return <- max(muhat.vals)
e.port.max<- efficient.portfolio(muhat.vals, cov.mat, target.return)
e.port.max

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = target.return)
## 
## Portfolio expected return:     0.00854 
## Portfolio standard deviation:  0.0198 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.9401 -0.4521 -0.1249  0.4215  0.3036 -0.0882

summary(e.port.max, risk.free=rf)

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = target.return)
## 
## Portfolio expected return:     0.00854 
## Portfolio standard deviation:  0.0198 
## Portfolio Sharpe Ratio:        0.41 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.9401 -0.4521 -0.1249  0.4215  0.3036 -0.0882

plot(e.port.max)

# Allowing shortsales...
plot(ef, plot.assets=TRUE, col="blue", lwd=2)
points(gmin.port$sd, gmin.port$er, col="orange", lwd=2)
text(0.005,0.005, labels="Glblminp (orange)", lwd=2)
points(e.port.max$sd, e.port.max$er, col="red", lwd=2)
text(e.port.max$sd, e.port.max$er, labels="MaxReturnp (red)", pos=4)
sr.ef = (e.port.max$er - rf)/e.port.max$sd
abline(a=rf, b=sr.ef, col="green", lwd=2)
abline(v=0, h=0)
points(0, rf, col="green", lwd=2)
text(0, rf, labels="rf", pos=4)
# Using the TBill+Tan Portfolio in Excel to get a line...
abline(a=rf, b=0.4321637, col="magenta", lwd=2)
points(0.00466, 0.0024256, col="cyan", lwd=2)
text(0.012, 0.0023, labels="TanTBillp (cyan)")
legend(x=0, y=0.017, c("Efficient Frontier", "Efficient Portfolio", "Tangent+T-Bills Portfolio"), col=c("blue", "green", "magenta"), lwd=2)

Tangency Portfolio [table in excel]

# Allowing short sales...
tan.port <- tangency.portfolio(muhat.vals, cov.mat, rf)
tan.port

## Call:
## tangency.portfolio(er = muhat.vals, cov.mat = cov.mat, risk.free = rf)
## 
## Portfolio expected return:     0.00288 
## Portfolio standard deviation:  0.0057 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.2800 -0.1172 -0.0436  0.0139  0.8881 -0.0212

summary(tan.port, risk.free=rf)

## Call:
## tangency.portfolio(er = muhat.vals, cov.mat = cov.mat, risk.free = rf)
## 
## Portfolio expected return:     0.00288 
## Portfolio standard deviation:  0.0057 
## Portfolio Sharpe Ratio:        0.433 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.2800 -0.1172 -0.0436  0.0139  0.8881 -0.0212

plot(tan.port)

VEURX, VEIEX, and VPACX have negative weights.

Tangency Portfolio + T-Bills (see magenta line above)

Annualized Mean, SD, and Sharpe Ratio for Tangency Portfolio [table in excel]

The annualized expected return, standard deviation and sharpe ratio values for the tangency portfolio are all larger than those of the global minimum variance portfolio, and the efficient portfolio with a target return of the maximum possible return on an asset.

Efficient Portfolio Frontier w/ No Short Sales

# No-short sales efficient frontiers...
gmin.port.ns = globalMin.portfolio(muhat.vals, cov.mat, shorts=FALSE)
tan.port.ns = tangency.portfolio(muhat.vals, cov.mat, rf, shorts=FALSE)
ef.ns = efficient.frontier(muhat.vals, cov.mat, alpha.min=0,                          
                        alpha.max=1, nport=20, shorts=FALSE)
plot(ef.ns, plot.assets=TRUE, col="blue", pch=16)
points(gmin.port.ns$sd, gmin.port.ns$er, col="orange")
text(tan.port.ns$sd, tan.port.ns$er,labels = "Tangent", pos = 2)
sr.tan.ns = (tan.port.ns$er - rf)/tan.port.ns$sd
text(0, rf, labels=expression(r[f]), pos = 4)
abline(a=rf, b=sr.tan.ns, col="green", lwd=2)
abline(h=0, v=0)
legend(x=0, y=0.0085, c("Efficient Frontier (NS)", "Global Min Portfolio (NS)", "Tangent Portfolio (NS)"), col=c("blue", "orange", "green"), lwd=2)

# Comparing no-short sales and short-sales frontiers...
# unrestricted efficient frontiers
plot(ef, plot.assets=TRUE, col="blue", pch=16)
sr.tan = (tan.port$er - rf)/tan.port$sd
text(0, rf, labels=expression(r[f]), pos = 3)
abline(a=rf, b=sr.tan, col="green", lwd=2)
abline(h=0, v=0)
legend(x=0, y=0.016, c("Eff Front", "Eff Front (NS)", "Tan Port", "Tan Port (NS)"), col=c("blue", "red", "green", "magenta"), lwd=2)
# short sales restricted frontiers
points(ef.ns$sd, ef.ns$er, type="b", pch=16, col="red")
abline(a=rf, b=sr.tan.ns, col="magenta", lwd=2)

Portfolio w/ Target Volatility = 0.02 [table in excel]

# Allowing Short Sales...
e.port.2<- efficient.portfolio(muhat.vals, cov.mat, 0.0086)
e.port.2

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = 0.0086)
## 
## Portfolio expected return:     0.0086 
## Portfolio standard deviation:  0.02 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.9466 -0.4554 -0.1257  0.4255  0.2979 -0.0889

summary(e.port.2, risk.free=rf)

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = 0.0086)
## 
## Portfolio expected return:     0.0086 
## Portfolio standard deviation:  0.02 
## Portfolio Sharpe Ratio:        0.409 
## Portfolio weights:
##   vfinx   veurx   veiex   vbltx   vbisx   vpacx 
##  0.9466 -0.4554 -0.1257  0.4255  0.2979 -0.0889

# No Short Sales...
e.port.2ns<- efficient.portfolio(muhat.vals, cov.mat, 0.0061, shorts=FALSE)
e.port.2ns

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = 0.0061, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.0061 
## Portfolio standard deviation:  0.02 
## Portfolio weights:
##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## 0.5088 0.0000 0.0000 0.4281 0.0631 0.0000

summary(e.port.2ns, risk.free=rf)

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = 0.0061, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.0061 
## Portfolio standard deviation:  0.02 
## Portfolio Sharpe Ratio:        0.284 
## Portfolio weights:
##  vfinx  veurx  veiex  vbltx  vbisx  vpacx 
## 0.5088 0.0000 0.0000 0.4281 0.0631 0.0000

For portfolios with a target volatilit of 2% per month, allowing for short sales results in an expected return equal to 0.86% per month. However, not allowing for short sales results in an expected return equal to 0.61% per month.

Tangency Portfolio with No Short Sales

tan.port.ns <- tangency.portfolio(muhat.vals, cov.mat, rf, shorts=FALSE)
summary(tan.port.ns, risk.free=rf)

## Call:
## tangency.portfolio(er = muhat.vals, cov.mat = cov.mat, risk.free = rf, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.00205 
## Portfolio standard deviation:  0.00527 
## Portfolio Sharpe Ratio:        0.31 
## Portfolio weights:
## vfinx veurx veiex vbltx vbisx vpacx 
## 0.145 0.000 0.000 0.000 0.855 0.000

plot(tan.port.ns)

The expected return, standard deviation, and sharpe ratio values for the tangency portfolio without shorts are smaller than those of the tangency portfolio with shorts.

Annualized Mean, SD, and Sharpe Ratio for Tangency Portfolio without Short Sales [table in excel]

The annualized expected return, mean, and Sharpe Ratio are all smaller than those of the tangency portfolio that allows for short sales.

ASSET ALLOCATION

Target Expexted Return 0.005 & No Shorts [table in excel]

e.port.6<- efficient.portfolio(er=muhat.vals, cov.mat, 0.005, shorts=FALSE)
e.port.6

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = 0.005, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.005 
## Portfolio standard deviation:  0.0159 
## Portfolio weights:
## vfinx veurx veiex vbltx vbisx vpacx 
## 0.411 0.000 0.000 0.310 0.279 0.000

summary(e.port.6, risk.free=rf)

## Call:
## efficient.portfolio(er = muhat.vals, cov.mat = cov.mat, target.return = 0.005, 
##     shorts = FALSE)
## 
## Portfolio expected return:     0.005 
## Portfolio standard deviation:  0.0159 
## Portfolio Sharpe Ratio:        0.288 
## Portfolio weights:
## vfinx veurx veiex vbltx vbisx vpacx 
## 0.411 0.000 0.000 0.310 0.279 0.000

plot(e.port.6)

About 41% of wealth is invested in VFINX , 31% in VBLTX, and 28% in VBISX for a target return of 0.5% monthly.

Portfolio w/ Target Return of 0.01 & No Shorts

#e.port.12<- efficient.portfolio(er=muhat.vals, cov.mat, 0.01, shorts=FALSE)
#e.port.12
#summary(e.port.12, risk.free=rf)
#plot(e.port.12)

There is no way to get a portfolio with a target return of 0.01 per month comprised of only risky assets and not allowing for short sales. The maximum expexted return on a portfolio comprised of only the risky assets would be about 0.0085 and all of the funds would be invested in VFINX.

REFERENCES

# https://en.wikipedia.org/wiki/2015%E2%80%9316_stock_market_selloff