library(tidyverse)
## -- Attaching packages -------------------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.2.1     v purrr   0.3.3
## v tibble  2.1.3     v dplyr   0.8.3
## v tidyr   1.0.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.4.0
## -- Conflicts ----------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo
## 
## Attaching package: 'xts'
## The following objects are masked from 'package:dplyr':
## 
##     first, last
library(PerformanceAnalytics)
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
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
library(quantmod)
## Loading required package: TTR
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(fPortfolio)
## Warning: package 'fPortfolio' was built under R version 3.6.3
## Loading required package: timeDate
## 
## Attaching package: 'timeDate'
## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness
## Loading required package: timeSeries
## Warning: package 'timeSeries' was built under R version 3.6.3
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
## Warning: package 'fBasics' was built under R version 3.6.3
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
## 
##     volatility
## Loading required package: fAssets
## Warning: package 'fAssets' was built under R version 3.6.3
library(readr)
library(fit.models)
## Warning: package 'fit.models' was built under R version 3.6.3
library(robust)
## Warning: package 'robust' was built under R version 3.6.3
library(Ecdat)
## Warning: package 'Ecdat' was built under R version 3.6.3
## Loading required package: Ecfun
## Warning: package 'Ecfun' was built under R version 3.6.3
## 
## Attaching package: 'Ecfun'
## The following object is masked from 'package:base':
## 
##     sign
## 
## Attaching package: 'Ecdat'
## The following object is masked from 'package:datasets':
## 
##     Orange
tickers <- c("SPY", "QQQ", "EEM","IWM","EFA","TLT","IYR","GLD")
data <- new.env()
getSymbols(tickers, src = 'yahoo', from = "2010-01-01", to = "2018-12-31", env = data, auto.assign = 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.
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## pausing 1 second between requests for more than 5 symbols
## [1] "SPY" "QQQ" "EEM" "IWM" "EFA" "TLT" "IYR" "GLD"

#Question 1

portfolioPrices <- NULL
for(ticker in tickers) {
  portfolioPrices <- cbind(portfolioPrices,
                           getSymbols(ticker, src = "yahoo", from='2010-01-01',to = "2018-12-31", periodicity = 'daily', auto.assign=FALSE)[,4])
}  
portfolioPrices <- portfolioPrices[apply(portfolioPrices,1,function(x) all(!is.na(x))),]
colnames(portfolioPrices) <- tickers
portfolioReturns <- na.omit(ROC(portfolioPrices, type="discrete"))
portfolioReturns <- as.timeSeries(portfolioReturns)
colnames(portfolioReturns) <- tickers
portfolioReturns <- as.timeSeries(portfolioReturns)
head(portfolioReturns)
## GMT
##                     SPY           QQQ          EEM           IWM
## 2010-01-05  0.002647093  0.0000000000  0.007258277 -0.0034386058
## 2010-01-06  0.000704057 -0.0060318615  0.002092073 -0.0009409818
## 2010-01-07  0.004221291  0.0006501734 -0.005799118  0.0073782887
## 2010-01-08  0.003327769  0.0082304747  0.007932804  0.0054542467
## 2010-01-11  0.001396552 -0.0040815898 -0.002083333 -0.0040296809
## 2010-01-12 -0.009326235 -0.0125108280 -0.016005636 -0.0108932614
##                      EFA           TLT          IYR           GLD
## 2010-01-05  0.0008813503  0.0064580894  0.002401157 -0.0009108014
## 2010-01-06  0.0042268581 -0.0133864256 -0.000435453  0.0164995902
## 2010-01-07 -0.0038583129  0.0016820139  0.008932440 -0.0061878037
## 2010-01-08  0.0079225530 -0.0004477891 -0.006694040  0.0049630301
## 2010-01-11  0.0082096245 -0.0054877141  0.004782630  0.0132889913
## 2010-01-12 -0.0117810116  0.0171170828 -0.016875854 -0.0209127164

#Question 2

or_matrix <- cor(portfolioReturns)
cov_matrix <- cov(portfolioReturns)
write.csv(cov_matrix, "covmatrix.csv")
mvp <- minvariancePortfolio(portfolioReturns, spec = portfolioSpec(), constraints = "Longonly")
mvpweights <- getWeights(mvp)
mvpweights
##       SPY       QQQ       EEM       IWM       EFA       TLT       IYR 
## 0.4259534 0.0000000 0.0000000 0.0000000 0.0000000 0.4495901 0.0000000 
##       GLD 
## 0.1244564

#Question 3

portfolioPrices <- NULL
for(ticker in tickers) {
  portfolioPrices <- cbind(portfolioPrices,
                           getSymbols.yahoo(ticker, from='2010-01-01',to = "2013-12-31", periodicity = 'weekly', auto.assign=FALSE)[,4])
}
portfolioPrices <- portfolioPrices[apply(portfolioPrices,1,function(x) all(!is.na(x))),]
colnames(portfolioPrices) <- tickers
portfolioReturns <- na.omit(ROC(portfolioPrices, type="discrete"))
portfolioReturns <- as.timeSeries(portfolioReturns)
colnames(portfolioReturns) <- tickers
portfolioReturns <- as.timeSeries(portfolioReturns)
head(portfolioReturns)
## GMT
##                    SPY          QQQ          EEM          IWM         EFA
## 2010-01-08  0.00648041  0.004765021 -0.009799393  0.007012701  0.02042254
## 2010-01-15 -0.02810409 -0.019400669 -0.045475991 -0.027236179 -0.05210490
## 2010-01-22 -0.02802146 -0.042646800 -0.044680254 -0.033089436 -0.03076081
## 2010-01-29 -0.01961866 -0.021354765 -0.027907028 -0.028627806 -0.03154930
## 2010-02-05  0.01587744  0.024636298  0.034821931  0.026084027  0.01221642
## 2010-02-12  0.02570986  0.027020839  0.020292859  0.039617001  0.02049806
##                     TLT          IYR          GLD
## 2010-01-08  0.013321336 -0.003886849  0.010918598
## 2010-01-15  0.018117621 -0.028181205 -0.041595966
## 2010-01-22 -0.007052973 -0.024760226 -0.008289094
## 2010-01-29  0.002950464 -0.021729208 -0.019815927
## 2010-02-05 -0.021355404  0.001870423  0.026444322
## 2010-02-12 -0.009463360  0.059976666  0.026603249
or_matrix <- cor(portfolioReturns)
cov_matrix <- cov(portfolioReturns)
write.csv(cov_matrix, "covmatrix.csv")
mvp <- minvariancePortfolio(portfolioReturns, spec = portfolioSpec(), constraints = "Longonly")
mvpweights <- getWeights(mvp)
mvpweights
##        SPY        QQQ        EEM        IWM        EFA        TLT 
## 0.49303564 0.00000000 0.00000000 0.00000000 0.00000000 0.48111726 
##        IYR        GLD 
## 0.00000000 0.02584709

#Question 4

portfolioPrices <- NULL
for(ticker in tickers) {
  portfolioPrices <- cbind(portfolioPrices,
                           getSymbols.yahoo(ticker, from='2010-01-01',to = "2013-12-31", periodicity = 'monthly', auto.assign=FALSE)[,4])
}
portfolioPrices <- portfolioPrices[apply(portfolioPrices,1,function(x) all(!is.na(x))),]
colnames(portfolioPrices) <- tickers
portfolioReturns <- na.omit(ROC(portfolioPrices, type="discrete"))
portfolioReturns <- as.timeSeries(portfolioReturns)
colnames(portfolioReturns) <- tickers
portfolioReturns <- as.timeSeries(portfolioReturns)
head(portfolioReturns)
## GMT
##                    SPY         QQQ          EEM         IWM          EFA
## 2010-02-01  0.03119470  0.04603872  0.017763846  0.04475126  0.002667664
## 2010-03-01  0.05652883  0.07596073  0.081108832  0.07961790  0.063854068
## 2010-04-01  0.01547007  0.02242529 -0.001661918  0.05678464 -0.028045731
## 2010-05-01 -0.07945455 -0.07392372 -0.093935817 -0.07536639 -0.111927954
## 2010-06-01 -0.05623116 -0.06337717 -0.020472390 -0.07743398 -0.037458651
## 2010-07-01  0.06830068  0.07258258  0.109324812  0.06380887  0.116104112
##                    TLT         IYR          GLD
## 2010-02-01 -0.00693316  0.05457055  0.032748217
## 2010-03-01 -0.02367185  0.08689957 -0.004386393
## 2010-04-01  0.02938544  0.06388108  0.058834366
## 2010-05-01  0.04743301 -0.05683531  0.030513141
## 2010-06-01  0.05440415 -0.05485489  0.023553189
## 2010-07-01 -0.01248154  0.09404794 -0.050871154
or_matrix <- cor(portfolioReturns)
cov_matrix <- cov(portfolioReturns)
write.csv(cov_matrix, "covmatrix.csv")
mvp <- minvariancePortfolio(portfolioReturns, spec = portfolioSpec(), constraints = "Longonly")
mvpweights <- getWeights(mvp)
mvpweights
##        SPY        QQQ        EEM        IWM        EFA        TLT 
## 0.44235680 0.00000000 0.00000000 0.04949215 0.00000000 0.48955621 
##        IYR        GLD 
## 0.00000000 0.01859485
retdata = read.csv("FFResearchDataFactors.CSV")
attach(retdata)
stocks5=cbind(N=Mkt.RF,SMB,HML)
fit = lm(cbind(Mkt.RF,SMB,HML)~RF)
options(digits=3)
head(SMB)
## [1] -2.3  -1.4  -1.32 0.04  -0.2  -0.04
## 758 Levels:  -0.01 -0.03 -0.04 -0.05 -0.07 -0.08 -0.09 -0.1 -0.11 ... SMB
pairs(cbind(Mkt.RF,SMB,HML))

cor(fit$residuals)
##         Mkt.RF    SMB     HML
## Mkt.RF 1.00000 0.1982 0.00895
## SMB    0.19824 1.0000 0.08089
## HML    0.00895 0.0809 1.00000
covRob(fit$residuals,cor=F)
## Call:
## covRob(data = fit$residuals, corr = F)
## 
## Robust Estimate of Covariance: 
##        Mkt.RF   SMB   HML
## Mkt.RF  66724 15796 -1920
## SMB     15796 48106  2864
## HML     -1920  2864 44271
## 
## Robust Estimate of Location: 
##  Mkt.RF     SMB     HML 
## -0.0452 -3.4413 -2.6780
cor.test(fit$residuals[,1], fit$residuals[,2])
## 
##  Pearson's product-moment correlation
## 
## data:  fit$residuals[, 1] and fit$residuals[, 2]
## t = 7, df = 1224, p-value = 2e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.144 0.251
## sample estimates:
##   cor 
## 0.198
cor.test(fit$residuals[,1], fit$residuals[,3])
## 
##  Pearson's product-moment correlation
## 
## data:  fit$residuals[, 1] and fit$residuals[, 3]
## t = 0.3, df = 1224, p-value = 0.8
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.0471  0.0649
## sample estimates:
##     cor 
## 0.00895
cor.test(fit$residuals[,2], fit$residuals[,3])
## 
##  Pearson's product-moment correlation
## 
## data:  fit$residuals[, 2] and fit$residuals[, 3]
## t = 3, df = 1224, p-value = 0.005
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.025 0.136
## sample estimates:
##    cor 
## 0.0809
pairs(fit$residuals)

n=dim(retdata)[1]
sigF = as.matrix(var(cbind(Mkt.RF,SMB,HML)))
sigF
##        Mkt.RF   SMB   HML
## Mkt.RF  65970 12827  1158
## SMB     12827 47529  4790
## HML      1158  4790 46390
bbeta = as.matrix(fit$coef)
bbeta = t( bbeta[-1,])
bbeta
##        RF-0.01 RF-0.02 RF-0.03 RF-0.06 RF0 RF0.01 RF0.02 RF0.03 RF0.04
## Mkt.RF     333     417     135      95 435    505    448    496    507
## SMB        276     545     218     207 339    414    422    351    320
## HML        407     212     214     105 229    374    376    376    338
##        RF0.05 RF0.06 RF0.07 RF0.08 RF0.09 RF0.1 RF0.11 RF0.12 RF0.13
## Mkt.RF    513    365    484    513    498   538    500    321    489
## SMB       475    280    424    237    339   320    331    346    321
## HML       401    335    363    347    262   341    311    209    344
##        RF0.14 RF0.15 RF0.16 RF0.17 RF0.18 RF0.19 RF0.2 RF0.21 RF0.22
## Mkt.RF    403    489    367    444    489    422   478    429    472
## SMB       386    356    270    365    346    255   265    299    345
## HML       271    355    298    289    235    347   278    268    304
##        RF0.23 RF0.24 RF0.25 RF0.26 RF0.27 RF0.28 RF0.29 RF0.3 RF0.31
## Mkt.RF    440    361    433    366    336    488    468   509    428
## SMB       401    291    279    272    310    336    249   314    306
## HML       354    324    341    260    361    387    268   386    316
##        RF0.32 RF0.33 RF0.34 RF0.35 RF0.36 RF0.37 RF0.38 RF0.39 RF0.4
## Mkt.RF    489    491    278    333    340    363    362    488   395
## SMB       330    502    356    447    397    353    304    326   342
## HML       314    339    336    398    398    271    332    280   358
##        RF0.41 RF0.42 RF0.43 RF0.44 RF0.45 RF0.46 RF0.47 RF0.48 RF0.49
## Mkt.RF    496    433    541    505    510    452    501    332    523
## SMB       402    331    343    334    428    337    311    253    413
## HML       265    314    313    259    417    353    397    443    266
##        RF0.5 RF0.51 RF0.52 RF0.53 RF0.54 RF0.55 RF0.56 RF0.57 RF0.58
## Mkt.RF   540    260    482    452    537    445    455    368    329
## SMB      375    299    286    297    442    202    432    254    372
## HML      493    402    368    403    292    440    321    431    433
##        RF0.59 RF0.6 RF0.61 RF0.62 RF0.63 RF0.64 RF0.65 RF0.66 RF0.67
## Mkt.RF    369   423    508    341    381    478    475    446    569
## SMB       462   283    378    302    590    422    240    223    233
## HML       397   388    314    376    338    273    289    274    278
##        RF0.68 RF0.69 RF0.7 RF0.71 RF0.72 RF0.73 RF0.74 RF0.75 RF0.76
## Mkt.RF    366    438   471    396     49    209    361    404    261
## SMB       331    365   393    209    372    211    444    191    259
## HML       285    123   322    520    635    489    551    278    597
##        RF0.77 RF0.78 RF0.79 RF0.8 RF0.81 RF0.82 RF0.83 RF0.86 RF0.87
## Mkt.RF    654    370    669   410    481    169    278    408    287
## SMB       554    414      3   367    412    277     22    277    356
## HML       389    276     73   425    445    418    113    644    279
##        RF0.89 RF0.92 RF0.95 RF0.96 RF0.98 RF0.99 RF1 RF1.02 RF1.04 RF1.05
## Mkt.RF     63    302    470    487    101    789  44    708    280    211
## SMB       144    404    621    329     16    591  98    614    615    433
## HML       397    710    186    436    636    242 387    692    721    356
##        RF1.06 RF1.07 RF1.08 RF1.1 RF1.13 RF1.15 RF1.2 RF1.21 RF1.24 RF1.26
## Mkt.RF    240    453    150   550    664    379   577    526    208    706
## SMB       402    250    684   661    489    536   569    509    196    446
## HML       493    366    552   291    218     37   667    257    372    434
##        RF1.28 RF1.31 RF1.35 RF1.49 RF1.54 RF1.57 RF1.6 RF1.65 RF1.66
## Mkt.RF    332    264    162    551    774    627   244    199    529
## SMB       154     25     72    312    529    338   626    691    337
## HML       675    576    684    312    517    703   412    729    617
##        RF1.81 RF1.82 RF10.38 RF10.54 RF11.24 RF14.71 RF2.13 RF2.14 RF2.41
## Mkt.RF    328     52     536     517     626     143    630    635    241
## SMB       261     95     607     751     714     739    487    330    318
## HML       336    325     184     518     230     597    686    165    167
##        RF2.46 RF2.66 RF2.73 RF2.9 RF2.95 RF2.98 RF3.12 RF3.14 RF3.51
## Mkt.RF    815     78    126   879    894    649    594    128    831
## SMB        96    219    349   718    712    152    263    216    743
## HML       133    297    742   529    497    738    392    318    595
##        RF3.54 RF3.56 RF3.83 RF3.84 RF3.9 RF3.93 RF4.21 RF4.39 RF4.66
## Mkt.RF    533    714    136    535   250    515    629    528    444
## SMB        82    683    535    159    88    609    724    715    344
## HML       750    310    531    487    61    726    332    163    173
##        RF4.68 RF4.75 RF4.76 RF4.8 RF4.86 RF5.08 RF5.12 RF5.21 RF5.26
## Mkt.RF    621    144    132   516    552    625    353    713    632
## SMB       530    281    580   380    234    528    610    439    339
## HML       270    515     60   520    335    596    730    619    526
##        RF5.47 RF5.6 RF5.8 RF5.89 RF6.16 RF6.35 RF6.52 RF6.58 RF6.93 RF7.18
## Mkt.RF    234   673   711    141    510    525    321    140    202    443
## SMB       157   442   531    122    356    721    158    161    233    527
## HML       136   430   744    648    745    521    593    340    527    375
##        RF7.72 RF7.81 RF8 RF8.37 RF8.8 RF9.85 RFRF
## Mkt.RF    631    134 243    620   540    311  895
## SMB       370    163  52    160   523    351  757
## HML       455    339 746    285   591    532  751
resig2 = apply((fit$resid)^2, 2, sum)/(n-3-1)
resig2 = diag(resig2)
cov_ff3 = sigF
cov_ff3
##        Mkt.RF   SMB   HML
## Mkt.RF  65970 12827  1158
## SMB     12827 47529  4790
## HML      1158  4790 46390
cov2cor(cov_ff3)
##        Mkt.RF   SMB    HML
## Mkt.RF 1.0000 0.229 0.0209
## SMB    0.2291 1.000 0.1020
## HML    0.0209 0.102 1.0000
cov_hist = cov(stocks5)
cov2cor(cov_hist)
##          N   SMB    HML
## N   1.0000 0.229 0.0209
## SMB 0.2291 1.000 0.1020
## HML 0.0209 0.102 1.0000
one.vec=rep(1,3)
a = solve(cov_ff3)%*%one.vec
b = t(one.vec)%*%a
mvp.w =a / as.numeric(b)
mvp.w
##         [,1]
## Mkt.RF 0.248
## SMB    0.337
## HML    0.415
a1 = solve(cov_hist)%*%one.vec
b1 = t(one.vec)%*%a1
mvp.w1 =a1 / as.numeric(b1)
mvp.w1
##      [,1]
## N   0.248
## SMB 0.337
## HML 0.415