rm(list = ls())
library(tidyquant)
## Loading required package: lubridate
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## Loading required package: PerformanceAnalytics
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
## Loading required package: quantmod
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## ══ Need to Learn tidyquant? ══════════════════════════════════════════════════════════════════════════════════════════════════════════════════
## Business Science offers a 1-hour course - Learning Lab #9: Performance Analysis & Portfolio Optimization with tidyquant!
## </> Learn more at: https://university.business-science.io/p/learning-labs-pro </>
library(lubridate)
library(timetk)
library(purrr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(quantmod)
library(quadprog)
library(lpSolve)
library(plyr)
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following object is masked from 'package:purrr':
## 
##     compact
  1. Import data
tickers<-c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "VNQ", "GLD")
data.env<-new.env()
l_ply(tickers, function(sym) try(getSymbols(sym, from="2010-01-01", env=data.env), silent = T))
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
## 
## This message is shown once per session and may be disabled by setting 
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
stocks<- tickers[tickers %in% ls(data.env)]
data<- xts()
for(i in seq_along(stocks)){
  symbol<-stocks[i]
  data <- merge(data, Ad(get(symbol, envir=data.env)))
}
data<- xts(coredata(data), order.by = as.Date(index(data)))
head(data)
##            SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## 2010-01-04     90.79006     41.51847     34.10928     54.61730     41.03314
## 2010-01-05     91.03041     41.51847     34.35684     54.42949     41.06930
## 2010-01-06     91.09449     41.26806     34.42872     54.37827     41.24289
## 2010-01-07     91.47905     41.29489     34.22906     54.77949     41.08377
## 2010-01-08     91.78343     41.63475     34.50059     55.07827     41.40926
## 2010-01-11     91.91164     41.46482     34.42872     54.85632     41.74921
##            TLT.Adjusted VNQ.Adjusted GLD.Adjusted
## 2010-01-04     66.13428     28.49396       109.80
## 2010-01-05     66.56142     28.46199       109.70
## 2010-01-06     65.67038     28.41083       111.51
## 2010-01-07     65.78078     28.71783       110.82
## 2010-01-08     65.75137     28.50676       111.37
## 2010-01-11     65.39054     28.67306       112.85
  1. Calculate weekly and monthly returns using log returns.
data.weekly<-data[endpoints(data, on="weeks", k= 1),]
lrets <- diff(log(data.weekly), lag= 5)
logretw <- na.omit(lrets)
head(logretw)
##            SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## 2010-02-12  -0.05868454  -0.06180639  -0.11674186 -0.055773408  -0.10100347
## 2010-02-19  -0.02224493  -0.02249731  -0.06322051 -0.009783898  -0.07476896
## 2010-02-26   0.01391211   0.01349534  -0.01654619  0.017185257  -0.02347730
## 2010-03-05   0.06192227   0.08185674   0.06742443  0.102828452   0.04653746
## 2010-03-12   0.07927751   0.09704303   0.10624681  0.133278216   0.07806817
## 2010-03-19   0.07493999   0.08286599   0.06909677  0.099591531   0.06761389
##             TLT.Adjusted VNQ.Adjusted GLD.Adjusted
## 2010-02-12  1.355619e-02 -0.059391410  -0.03965541
## 2010-02-19 -1.453119e-02 -0.007428305  -0.01261760
## 2010-02-26 -6.745915e-05  0.046551885   0.02086874
## 2010-03-05 -1.565164e-02  0.093653053   0.04475535
## 2010-03-12 -1.327378e-02  0.133056202   0.03076005
## 2010-03-19  1.496284e-02  0.155170421   0.01151785
data.monthly<-data[endpoints(data, on="months", k= 1),]

lrets <- diff(log(data.monthly), lag=5)
logretm <- na.omit(lrets)
head(logretm)
##            SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## 2010-06-30 -0.030748107  0.003053435  -0.01879835   0.01915032  -0.10342144
## 2010-07-30  0.004603132  0.028112328   0.06734515   0.04051210   0.00375850
## 2010-08-31 -0.100519229 -0.098830590  -0.04354459  -0.11593720  -0.09682894
## 2010-09-30 -0.030101864  0.002738416   0.06927834  -0.05387136   0.02667301
## 2010-10-29  0.090177343  0.140988635   0.19763202   0.06513086   0.18272455
## 2010-11-30  0.143304749  0.200879068   0.18223220   0.17998396   0.15411990
##            TLT.Adjusted VNQ.Adjusted GLD.Adjusted
## 2010-06-30   0.11464143   0.11204475   0.13833299
## 2010-07-30   0.10856383   0.14932334   0.05389885
## 2010-08-31   0.20992471   0.03947716   0.11378755
## 2010-09-30   0.15171975   0.01415243   0.10326924
## 2010-10-29   0.05617549   0.11531152   0.10937330
## 2010-11-30  -0.01722054   0.15015252   0.10698640
  1. Convert monthly returns into tibble format.
tibble::as_tibble(logretm, rownames = 'date')
## # A tibble: 131 x 9
##    date  SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
##    <chr>        <dbl>        <dbl>        <dbl>        <dbl>        <dbl>
##  1 2010…     -0.0307       0.00305      -0.0188       0.0192     -0.103  
##  2 2010…      0.00460      0.0281        0.0673       0.0405      0.00376
##  3 2010…     -0.101       -0.0988       -0.0435      -0.116      -0.0968 
##  4 2010…     -0.0301       0.00274       0.0693      -0.0539      0.0267 
##  5 2010…      0.0902       0.141         0.198        0.0651      0.183  
##  6 2010…      0.143        0.201         0.182        0.180       0.154  
##  7 2010…      0.142        0.177         0.149        0.192       0.124  
##  8 2011…      0.211        0.258         0.142        0.266       0.184  
##  9 2011…      0.159        0.165         0.0309       0.202       0.123  
## 10 2011…      0.122        0.0993        0.0622       0.187       0.0618 
## # … with 121 more rows, and 3 more variables: TLT.Adjusted <dbl>,
## #   VNQ.Adjusted <dbl>, GLD.Adjusted <dbl>
  1. Download Fama French 3 factors data and change to digit numbers (not in percetage).
library(readxl)
f.f <- read_excel("F-F_Research_Data_Factors.xlsx")
## New names:
## * `` -> ...1
F_F_Research_Data_Factors <- read_excel("F-F_Research_Data_Factors.xlsx")
## New names:
## * `` -> ...1
date <- seq(as.Date("1926-06-30"), length=1136, by="1 month")
f.f <- xts(coredata(F_F_Research_Data_Factors[, -1])/100, order.by = date)
names(f.f)[names(f.f) == "Mkt-RF"] <- "MKT"
head(f.f)
##                MKT     SMB     HML     RF
## 1926-06-30  0.0296 -0.0230 -0.0287 0.0022
## 1926-07-30  0.0264 -0.0140  0.0419 0.0025
## 1926-08-30  0.0036 -0.0132  0.0001 0.0023
## 1926-09-30 -0.0324  0.0004  0.0051 0.0032
## 1926-10-30  0.0253 -0.0020 -0.0035 0.0031
## 1926-11-30  0.0262 -0.0004 -0.0002 0.0028
  1. Merge monthly return data in question 3 and 4 into tibble format.
monthly.ret <- tk_tbl(logretm)
ff3 <- tk_tbl(f.f)
merged.monthly.ret <- monthly.ret %>%
  inner_join(ff3)
## Joining, by = "index"
head(merged.monthly.ret)
## # A tibble: 6 x 13
##   index      SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
##   <date>            <dbl>        <dbl>        <dbl>        <dbl>        <dbl>
## 1 2010-06-30     -0.0307       0.00305      -0.0188       0.0192     -0.103  
## 2 2010-07-30      0.00460      0.0281        0.0673       0.0405      0.00376
## 3 2010-09-30     -0.0301       0.00274       0.0693      -0.0539      0.0267 
## 4 2010-11-30      0.143        0.201         0.182        0.180       0.154  
## 5 2011-06-30      0.0345       0.0222        0.0483       0.0624      0.0310 
## 6 2011-09-30     -0.177       -0.114        -0.344       -0.288      -0.265  
## # … with 7 more variables: TLT.Adjusted <dbl>, VNQ.Adjusted <dbl>,
## #   GLD.Adjusted <dbl>, MKT <dbl>, SMB <dbl>, HML <dbl>, RF <dbl>
  1. Based on CAPM model, compute MVP monthly returns based on estimated covariance matrix for the 8-asset portfolio by using past 60-month returns from 2015/01 - 2021/03.
ret4 <- merged.monthly.ret %>% 
  select(-RF) %>% 
  gather(key = stock, value = price, -index) %>%
  group_by(stock) %>% 
  tq_transmute( mutate_fun = periodReturn, 
                period     = "monthly", 
                type       = "arithmetic",
                col_rename = "monthly.returns") %>% 
  ungroup() %>% 
  spread(stock, monthly.returns) %>% 
  bind_cols(., RF = merged.monthly.ret$RF) %>%
  mutate(RF = RF/12) %>%
  filter(index >= as.Date("2014-09-30") & index <= as.Date("2021-03-01")) %>%
  mutate(SPY_RF = SPY.Adjusted - RF,
         QQQ_RF = QQQ.Adjusted - RF,
         EEM_RF = EEM.Adjusted - RF,
         IWM_RF = IWM.Adjusted - RF,
         EFA_RF = EFA.Adjusted - RF,
         TLT_RF = TLT.Adjusted - RF,
         VNQ_RF = VNQ.Adjusted - RF,
         GLD_RF = GLD.Adjusted - RF,
         MKT_RF = MKT) %>%
  slice(-1) %>%
  select(index, ends_with("_RF"))

head(ret4)
## # A tibble: 6 x 10
##   index      SPY_RF QQQ_RF   EEM_RF  IWM_RF  EFA_RF TLT_RF  VNQ_RF GLD_RF MKT_RF
##   <date>      <dbl>  <dbl>    <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>  <dbl>
## 1 2015-01-30 -0.929 -0.849 -10.7     -0.945  0.943   1.78   38.6   -0.950  1.44 
## 2 2015-04-30  3.14   0.101  -1.36   -48.3   -1.64   -0.763  -0.966 -4.64  -0.779
## 3 2015-06-30  1.55   1.93   -0.518    0.780  0.0991 -4.99  -29.5   -8.80   0.132
## 4 2015-09-30 -2.81  -1.83  -12.5     -2.22  -3.48   -0.944  -0.775 -0.362  4.03 
## 5 2015-10-30 -0.929 -1.66   -0.397   -0.336 -0.470  -2.32   -2.00  -0.284 -0.928
## 6 2015-11-30 -4.58   0.916  -0.0266  -0.377 -0.345   2.90    1.41   1.27  -4.88
ret4.reg <- ret4 %>% 
  lm(formula = cbind(SPY_RF, QQQ_RF, EEM_RF, IWM_RF, EFA_RF, TLT_RF, VNQ_RF, GLD_RF) ~ MKT_RF, data =.)
head(ret4.reg)
## $coefficients
##                 SPY_RF   QQQ_RF     EEM_RF     IWM_RF   EFA_RF     TLT_RF
## (Intercept) -0.4850552 11.52168 -2.7229103 -2.8156904 5.448528 -0.9686981
## MKT_RF      -0.1380028  3.95935 -0.6318787  0.5795739 1.934063 -0.3102419
##                VNQ_RF      GLD_RF
## (Intercept) 0.6115007 -1.37465840
## MKT_RF      0.3857227  0.02083774
## 
## $residuals
##         SPY_RF      QQQ_RF      EEM_RF      IWM_RF      EFA_RF      TLT_RF
## 1  -0.24590276 -18.0586455 -7.02334589   1.0379367  -7.2834567  3.19675397
## 2   3.51739874  -8.3382404  0.86849778 -44.9868705  -5.5779060 -0.03545672
## 3   2.05570168 -10.1202703  2.28902844   3.5187187  -5.6053638 -3.98376437
## 4  -1.77166015 -29.3163757 -7.24779018  -1.7455090 -16.7302298  1.27576023
## 5  -0.57173222  -9.5123651  1.73956089   3.0174354  -4.1246141 -1.64370899
## 6  -4.76341736   8.6959089 -0.38409264   5.2641513   3.6348448  2.35179457
## 7   3.78385605  -0.7758745 -0.87279118   0.6240853  -1.5721338  1.29104244
## 8  -0.07511297  -3.3366853  1.51420413   3.6389052  -2.3062382  0.22546252
## 9   0.19393215  -4.2431608  0.87239812   4.2514587  -2.2378553 -2.53748820
## 10  0.21162611 -12.4174906  1.08032009   2.5775319  -6.2898741  1.17748257
## 11  1.23040173  -9.2356509 -3.11788579   2.2210424  14.9889309 -0.36243802
## 12  0.80462433  -9.6406913  2.45077514   4.1808101  -4.9217323  0.03734477
## 13 -0.23345656 -18.1299216  2.97376583   0.9884631  -9.0812302 -3.29856115
## 14  1.69365375   5.4192701 -2.81836869 -11.8295517  -0.8668929 -0.36723767
## 15  3.09318767 -21.5652026  0.01252498   1.3127008  -7.0367459 -2.01883700
## 16 -0.15248254  -7.3646968  0.25173326   1.9346224  -4.1772857  1.60720936
## 17  0.02043872 -14.2580011  3.11267489   3.0975418  -6.8153493  1.06194770
## 18 -3.90626240 171.4224261  4.78648279   6.3370151  71.9422060  1.65209221
## 19 -0.03739648 -10.8382832  2.10230816   2.3092368  -6.0072088  0.66752548
## 20 -5.38317961  -5.2199662 -3.30020405   1.2388789  -5.2939725 -0.57654472
## 21 -0.94898633  15.5536732 -1.74155481   6.5498137   7.2960163 -0.08024892
## 22  1.48476846  -8.7197563  2.45175871   4.4615830  -1.9339085  0.35986993
##          VNQ_RF     GLD_RF
## 1   37.45450197  0.3944690
## 2   -1.27761112 -3.2457296
## 3  -30.11789175 -7.4270328
## 4   -2.94175373  0.9288641
## 5   -2.24904352  1.1100783
## 6    2.67964040  2.7473262
## 7    1.74939850 -1.2840954
## 8   -0.56917002  0.5013714
## 9   -2.10813837 -7.3361966
## 10  -0.69511776  1.6883273
## 11  -1.91420422  0.2453766
## 12  -0.55881275  1.6016561
## 13  -5.45097649  1.4461135
## 14  -0.09583996 -0.4644368
## 15   2.20847893 -1.4203897
## 16  -0.23766522  3.8282189
## 17  -0.48102175  1.1953462
## 18  -4.31052179  1.3780657
## 19  -0.74111744  1.1804661
## 20  -1.42196312  1.3604021
## 21   1.18734503  1.2374916
## 22   9.89148418  0.3343078
## 
## $effects
##                 SPY_RF      QQQ_RF      EEM_RF      IWM_RF      EFA_RF
## (Intercept)  1.7176637 -38.0481352 10.21917816  15.5478825 -17.7434303
## MKT_RF      -1.5103888  43.3335865 -6.91567277   6.3432169  21.1675870
##              2.2853157  -7.2171276  3.63334955   0.9716749  -4.5504132
##             -0.3093854 -28.2012255 -5.17551454 -19.9298597 -17.1632029
##             -0.6771692  -6.1232259  2.88601580   4.7207869  -2.6652291
##             -6.1164151  13.8946535  0.02560844  22.7938290   6.6001392
##              3.0802692   3.4808900 -0.07957694   9.9154559   0.6092678
##             -0.3650171   0.3200262  2.55172096   7.6823729  -0.6241860
##             -0.2191073  -0.4078397  1.83719673   9.8569982  -0.4071686
##              0.4202480  -9.4838983  2.41224415   0.2967911  -5.2095841
##              1.4067805  -6.2552896 -1.80500306   0.3493316  16.1081409
##              0.8555053  -6.4782935  3.68954445   3.9011413  -3.6510360
##              0.4284127 -15.8537714  4.57335762  -7.0420906  -8.5480474
##              0.4263792  10.4936720 -2.35804333   4.6126578   1.9949265
##              3.1945925 -18.4760905  1.28113146   0.3920951  -5.8270361
##             -0.3458815  -3.8479675  1.34624175   4.7538432  -2.6117233
##              0.3474638 -11.4961544  4.51452261  -0.6852400  -5.8779818
##             -3.1784200 173.6028812  6.42503540  -2.5304622  72.3957537
##              0.1611482  -7.8900740  3.42828111   0.1563327  -4.9147549
##             -5.6264839  -1.6308484 -2.23516744   4.6911899  -3.6681702
##             -2.9540409  21.6982359 -1.71692948  32.3513646  11.0483974
##              1.4738431  -5.4677077  3.65402800   4.9659771  -0.5886069
##                  TLT_RF      VNQ_RF     GLD_RF
## (Intercept)  3.29040946  -1.3101086  6.5318911
## MKT_RF      -3.39548004   4.2215891  0.2280611
##             -4.57680868 -37.1113822 -7.6698997
##              0.47595694 -12.6562015 -0.4616519
##             -2.18055390  -8.5029472  1.1791553
##              2.02420788  -0.8204158  3.9779247
##              0.85452754  -3.1841579 -0.6581205
##             -0.28044100  -6.4158836  0.7421936
##             -3.02273776  -7.6830449 -6.9807312
##              0.58795936  -7.6422704  1.4650048
##             -0.94655298  -8.7901839  0.0520736
##             -0.52571995  -7.1577704  1.5251958
##             -3.96410932 -13.3986208  0.8008023
##             -0.70920304  -3.7851202  0.6863506
##             -2.59037629  -4.5020044 -1.5438895
##              1.08511865  -6.2974030  3.9791915
##              0.45256426  -7.6895358  0.8617863
##              0.97547811 -12.4037943  0.6713312
##              0.07969254  -7.6660259  0.9665258
##             -1.09026460  -7.3715405  1.5578383
##             -0.29846357  -0.8733704  3.0751773
##             -0.19282778   3.4289568  0.3153913
## 
## $rank
## [1] 2
## 
## $fitted.values
##         SPY_RF      QQQ_RF     EEM_RF     IWM_RF       EFA_RF      TLT_RF
## 1  -0.68329740  17.2093214 -3.6306090 -1.9831278  8.226824815 -1.41436306
## 2  -0.37761979   8.4393232 -2.2309917 -3.2668896  3.942856761 -0.72717429
## 3  -0.50332031  12.0457153 -2.8065413 -2.7389821  5.704506887 -1.00975954
## 4  -1.04154719  27.4876342 -5.2709405 -0.4785774 13.247574229 -2.21973849
## 5  -0.35702421   7.8484286 -2.1366900 -3.3533854  3.654216605 -0.68087369
## 6   0.18770861  -7.7801478  0.3574983 -5.6411132 -3.980029185  0.54373113
## 7  -0.09584907   0.3552218 -0.9408377 -4.4502491 -0.006064475 -0.09373018
## 8  -0.27647879   5.5375521 -1.7678937 -3.6916540  2.525399825 -0.49980087
## 9  -0.22271320   3.9949984 -1.5217152 -3.9174546  1.771893338 -0.37893134
## 10 -0.49415432  11.7827397 -2.7645726 -2.7774767  5.576048597 -0.98915362
## 11 -0.48007574  11.3788205 -2.7001106 -2.8366028  5.378742181 -0.95750382
## 12 -0.42527860   9.8066711 -2.4492088 -3.0667358  4.610778744 -0.83431526
## 13 -0.69205948  17.4607088 -3.6707283 -1.9463295  8.349622466 -1.43406096
## 14  0.15027859  -6.7062675  0.1861161 -5.4839174 -3.455460210  0.45958534
## 15 -0.44733927  10.4395995 -2.5502188 -2.9740870  4.919951578 -0.88390949
## 16 -0.31861665   6.7465021 -1.9608318 -3.5146866  3.115947651 -0.59453029
## 17 -0.54585368  13.2660127 -3.0012904 -2.5603536  6.300597675 -1.10537811
## 18 -0.72086590  18.2871750 -3.8026253 -1.8253505  8.753334646 -1.49882019
## 19 -0.48975425  11.6565003 -2.7444259 -2.7959558  5.514383152 -0.97926190
## 20 -0.29682605   6.1213224 -1.8610584 -3.6062011  2.810559876 -0.54554322
## 21  0.47242160 -15.9486640  1.6611241 -6.8368294 -7.970185358  1.18378971
## 22 -0.39829161   9.0324051 -2.3256425 -3.1800736  4.232565333 -0.77364627
##         VNQ_RF    GLD_RF
## 1   1.16559442 -1.344725
## 2   0.31121490 -1.390881
## 3   0.66255220 -1.371900
## 4   2.16691504 -1.290631
## 5   0.25364952 -1.393990
## 6  -1.26889762 -1.476242
## 7  -0.47634407 -1.433427
## 8   0.02852227 -1.406152
## 9  -0.12175442 -1.414271
## 10  0.63693293 -1.373284
## 11  0.59758284 -1.375410
## 12  0.44442291 -1.383684
## 13  1.19008476 -1.343402
## 14 -1.16427943 -1.470591
## 15  0.50608325 -1.380353
## 16  0.14629905 -1.399790
## 17  0.78143446 -1.365478
## 18  1.27059969 -1.339052
## 19  0.62463459 -1.373949
## 20  0.08539359 -1.403080
## 21 -2.06468035 -1.519233
## 22  0.36899336 -1.387759
## 
## $assign
## [1] 0 1
b_hat <- ret4.reg$coefficients
diagD_hat <- ret4.reg$residuals %>% cov() %>% diag() %>% diag(nrow = length(.))
diagD_hat
##          [,1]     [,2]     [,3]     [,4]     [,5]    [,6]     [,7]     [,8]
## [1,] 5.678079    0.000 0.000000   0.0000   0.0000 0.00000   0.0000 0.000000
## [2,] 0.000000 1562.607 0.000000   0.0000   0.0000 0.00000   0.0000 0.000000
## [3,] 0.000000    0.000 9.787619   0.0000   0.0000 0.00000   0.0000 0.000000
## [4,] 0.000000    0.000 0.000000 114.4623   0.0000 0.00000   0.0000 0.000000
## [5,] 0.000000    0.000 0.000000   0.0000 296.2076 0.00000   0.0000 0.000000
## [6,] 0.000000    0.000 0.000000   0.0000   0.0000 3.24153   0.0000 0.000000
## [7,] 0.000000    0.000 0.000000   0.0000   0.0000 0.00000 119.0491 0.000000
## [8,] 0.000000    0.000 0.000000   0.0000   0.0000 0.00000   0.0000 7.804049
sigm2 = var(ret4$MKT_RF)
Omega = sigm2*t(b_hat)%*%b_hat + diagD_hat
Omega
##            SPY_RF     QQQ_RF     EEM_RF      IWM_RF    EFA_RF     TLT_RF
## SPY_RF   7.128748  -34.99456    8.03107    7.334148 -16.59728   2.924381
## QQQ_RF -34.994565 2409.23241 -193.22042 -171.958159 401.75713 -70.669510
## EEM_RF   8.031070 -193.22042   54.35616   41.643173 -91.59509  16.163600
## IWM_RF   7.334148 -171.95816   41.64317  161.600587 -81.11386  14.532431
## EFA_RF -16.597277  401.75713  -91.59509  -81.113859 486.87670 -33.528352
## TLT_RF   2.924381  -70.66951   16.16360   14.532431 -33.52835   9.143073
## VNQ_RF  -1.995513   48.89914  -10.88781   -8.546023  23.25986  -4.061426
## GLD_RF   3.786963  -89.87206   21.27551   22.146996 -42.49257   7.558783
##            VNQ_RF     GLD_RF
## SPY_RF  -1.995513   3.786963
## QQQ_RF  48.899143 -89.872064
## EEM_RF -10.887813  21.275506
## IWM_RF  -8.546023  22.146996
## EFA_RF  23.259860 -42.492569
## TLT_RF  -4.061426   7.558783
## VNQ_RF 122.030675  -4.748991
## GLD_RF  -4.748991  18.585359
one.vec = rep(1, 8)
sigma.inv.mat = solve(Omega)
top.mat = sigma.inv.mat%*%one.vec
bot.val = as.numeric ((t(one.vec)%*%sigma.inv.mat%*%one.vec))
m.mat = top.mat/bot.val
m.mat[,1]
##       SPY_RF       QQQ_RF       EEM_RF       IWM_RF       EFA_RF       TLT_RF 
##  0.395033907  0.012161506 -0.086550510 -0.007174135  0.035555578  0.483413745 
##       VNQ_RF       GLD_RF 
##  0.031856219  0.135703689
  1. Based on FF 3-factor model, compute MVP monthly returns covariance matrix for the 8-asset portfolio by using past 60-month returns from 2015/01 - 2021/03.
N <- dim(merged.monthly.ret[,-13])[1]
fit3 = lm(formula = cbind(SPY.Adjusted,QQQ.Adjusted,EEM.Adjusted,IWM.Adjusted,EFA.Adjusted,TLT.Adjusted,VNQ.Adjusted,GLD.Adjusted) ~ MKT + SMB + HML, data = merged.monthly.ret)
sigF3 = as.matrix(var(cbind(merged.monthly.ret$MKT, 
                            merged.monthly.ret$SMB, 
                            merged.monthly.ret$HML)))
bbeta3 = as.matrix(fit3$coefficients)
bbeta3 = bbeta3[-1,]
bbeta3
##     SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## MKT   -0.7295811   -0.5119661   -1.0230772   -0.6575093   -0.8379082
## SMB    0.6443993    0.6839078    0.9904405    0.4400724    0.5174659
## HML    0.9168393    0.7802921    1.2032259    1.3718865    1.1658264
##     TLT.Adjusted VNQ.Adjusted GLD.Adjusted
## MKT   0.41756813   -0.9222459    0.1498878
## SMB  -0.05662353    0.3830940    0.1556907
## HML  -0.18012818    0.9853922    0.5694025
sigeps3 = crossprod(fit3$residuals)/(N-5)
sigeps3 = diag(diag(sigeps3))
cov_3f = t(bbeta3) %*% sigF3 %*% (bbeta3) + sigeps3
cov_3f
##               SPY.Adjusted  QQQ.Adjusted  EEM.Adjusted  IWM.Adjusted
## SPY.Adjusted  0.0050222846  0.0009706178  0.0017327336  0.0013983983
## QQQ.Adjusted  0.0009706178  0.0058528265  0.0013444753  0.0010813735
## EEM.Adjusted  0.0017327336  0.0013444753  0.0171420472  0.0018938223
## IWM.Adjusted  0.0013983983  0.0010813735  0.0018938223  0.0122029057
## EFA.Adjusted  0.0014871211  0.0011143221  0.0020259801  0.0017010682
## TLT.Adjusted -0.0005250839 -0.0003355795 -0.0007139824 -0.0005402463
## VNQ.Adjusted  0.0014735265  0.0010550075  0.0020040166  0.0016449981
## GLD.Adjusted  0.0001423906  0.0001953799  0.0001872847  0.0002967451
##               EFA.Adjusted  TLT.Adjusted  VNQ.Adjusted GLD.Adjusted
## SPY.Adjusted  0.0014871211 -0.0005250839  1.473527e-03 1.423906e-04
## QQQ.Adjusted  0.0011143221 -0.0003355795  1.055008e-03 1.953799e-04
## EEM.Adjusted  0.0020259801 -0.0007139824  2.004017e-03 1.872847e-04
## IWM.Adjusted  0.0017010682 -0.0005402463  1.644998e-03 2.967451e-04
## EFA.Adjusted  0.0097143400 -0.0006529375  1.800333e-03 1.598160e-04
## TLT.Adjusted -0.0006529375  0.0122836560 -7.310754e-04 1.047757e-04
## VNQ.Adjusted  0.0018003326 -0.0007310754  1.132153e-02 3.154319e-05
## GLD.Adjusted  0.0001598160  0.0001047757  3.154319e-05 9.597635e-03
sigma.inv.mat.ff3 = solve(cov_3f)
top.mat.ff3 = sigma.inv.mat.ff3%*%one.vec
bot.val.ff3 = as.numeric ((t(one.vec)%*%sigma.inv.mat.ff3%*%one.vec))
m.mat.ff3 = top.mat.ff3/bot.val.ff3
m.mat.ff3[,1]
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted TLT.Adjusted 
##   0.22049242   0.19548700   0.03600436   0.06645344   0.08668738   0.16064570 
## VNQ.Adjusted GLD.Adjusted 
##   0.07749151   0.15673820