#rm(list=ls())
library(quantmod)
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## 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(plyr)
library(fBasics)
## Loading required package: timeDate
## Loading required package: timeSeries
##
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
##
## time<-
##
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
##
## volatility
library(tidyverse)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
## Registered S3 method overwritten by 'rvest':
## method from
## read_xml.response xml2
## ── Attaching packages ────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.1 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.5
## ✔ tidyr 1.0.2 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::arrange() masks plyr::arrange()
## ✖ purrr::compact() masks plyr::compact()
## ✖ dplyr::count() masks plyr::count()
## ✖ dplyr::failwith() masks plyr::failwith()
## ✖ dplyr::filter() masks timeSeries::filter(), stats::filter()
## ✖ dplyr::first() masks xts::first()
## ✖ dplyr::id() masks plyr::id()
## ✖ dplyr::lag() masks timeSeries::lag(), stats::lag()
## ✖ dplyr::last() masks xts::last()
## ✖ dplyr::mutate() masks plyr::mutate()
## ✖ dplyr::rename() masks plyr::rename()
## ✖ dplyr::summarise() masks plyr::summarise()
## ✖ dplyr::summarize() masks plyr::summarize()
library(openair)
con = gzcon(url('https://github.com/systematicinvestor/SIT/raw/master/sit.gz', 'rb'))
source(con)
close(con)
library(xts)
#Download ETF daily data from yahoo with ticker names of SPY, QQQ, EEM, IWM, EFA, TLT, IYR, GLD from 2010 to 2018.
tickers<-c("SPY", "QQQ", "EEM", "IWM", "EFA", "TLT", "IYR", "GLD" )
data.env<-new.env()
getSymbols(tickers ,from = "2010/01/01", to = "2018/12/31", env = data.env, auto.assign = TRUE)
## '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"
l_ply(tickers, function(sym) try(getSymbols(sym ,from = "2010/01/01" ,to = "2018/12/31", env=data.env, auto.assign = TRUE), silent=T))
stocks <- tickers[tickers %in% ls(data.env)]
stocks
## [1] "SPY" "QQQ" "EEM" "IWM" "EFA" "TLT" "IYR" "GLD"
data <- xts()
for(i in seq_along(stocks))
{
symbol <- stocks[i]
data <- merge(data, Ad(get(symbol, envir=data.env)))}
head(data)
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2010-01-04 92.24605 41.71869 34.66187 55.28631
## 2010-01-05 92.49020 41.71869 34.91345 55.09619
## 2010-01-06 92.55533 41.46705 34.98648 55.04436
## 2010-01-07 92.94606 41.49401 34.78360 55.45049
## 2010-01-08 93.25535 41.83553 35.05952 55.75293
## 2010-01-11 93.38558 41.66477 34.98648 55.52828
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2010-01-04 42.01963 67.05941 31.17820 109.80
## 2010-01-05 42.05667 67.49248 31.25306 109.70
## 2010-01-06 42.23443 66.58901 31.23945 111.51
## 2010-01-07 42.07148 66.70103 31.51849 110.82
## 2010-01-08 42.40480 66.67115 31.30751 111.37
## 2010-01-11 42.75292 66.30531 31.45724 112.85
tail(data)
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2018-12-20 239.7584 150.2941 37.69136 129.5411
## 2018-12-21 234.8459 145.6359 37.42901 126.1795
## 2018-12-24 228.6406 142.0244 37.07921 123.7320
## 2018-12-26 240.1926 150.8923 37.81768 129.6788
## 2018-12-27 242.0367 151.4763 37.76910 130.2194
## 2018-12-28 241.7244 151.3971 38.12862 130.5929
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2018-12-20 56.70416 117.4275 72.57053 119.24
## 2018-12-21 55.71682 117.4080 71.66856 118.72
## 2018-12-24 55.06827 117.9916 69.08737 120.02
## 2018-12-26 56.34601 116.7272 71.35190 119.66
## 2018-12-27 56.51057 116.7467 71.46705 120.57
## 2018-12-28 56.81064 117.7290 71.68774 121.06
#Q1. Use the data to calculate daily returns. Find their means, median, 1-st quartile, third-quartile, standard deviation, skewness and kurtosis.
SPYreturn = dailyReturn(data$SPY.Adjusted)
QQQreturn = dailyReturn(data$QQQ.Adjusted)
EEMreturn = dailyReturn(data$EEM.Adjusted)
IWMreturn = dailyReturn(data$IWM.Adjusted)
EFAreturn = dailyReturn(data$EFA.Adjusted)
TLTreturn = dailyReturn(data$TLT.Adjusted)
IYRreturn = dailyReturn(data$IYR.Adjusted)
GLDreturn = dailyReturn(data$GLD.Adjusted)
dailyreturns = merge(SPYreturn,QQQreturn,EEMreturn,IWMreturn,EFAreturn,TLTreturn,IYRreturn,GLDreturn )
colnames(dailyreturns)
## [1] "daily.returns" "daily.returns.1" "daily.returns.2" "daily.returns.3"
## [5] "daily.returns.4" "daily.returns.5" "daily.returns.6" "daily.returns.7"
names(dailyreturns)[names(dailyreturns) =="daily.returns" ] <- "SPY"
names(dailyreturns)[names(dailyreturns) =="daily.returns.1" ] <- "QQQ"
names(dailyreturns)[names(dailyreturns) =="daily.returns.2" ] <- "EEM"
names(dailyreturns)[names(dailyreturns) =="daily.returns.3" ] <- "IWM"
names(dailyreturns)[names(dailyreturns) =="daily.returns.4" ] <- "EFA"
names(dailyreturns)[names(dailyreturns) =="daily.returns.5" ] <- "TLT"
names(dailyreturns)[names(dailyreturns) =="daily.returns.6" ] <- "IYR"
names(dailyreturns)[names(dailyreturns) =="daily.returns.7" ] <- "GLD"
head(dailyreturns)
## SPY QQQ EEM IWM
## 2010-01-04 0.0000000000 0.0000000000 0.000000000 0.000000000
## 2010-01-05 0.0026467909 0.0000000000 0.007258004 -0.003438772
## 2010-01-06 0.0007041178 -0.0060318290 0.002091916 -0.000940791
## 2010-01-07 0.0042216046 0.0006501065 -0.005798954 0.007378195
## 2010-01-08 0.0033276074 0.0082306344 0.007932619 0.005454398
## 2010-01-11 0.0013965312 -0.0040817220 -0.002083314 -0.004029546
## EFA TLT IYR GLD
## 2010-01-04 0.0000000000 0.0000000000 0.0000000000 0.0000000000
## 2010-01-05 0.0008813261 0.0064579602 0.0024009402 -0.0009108014
## 2010-01-06 0.0042267971 -0.0133861586 -0.0004354454 0.0164995902
## 2010-01-07 -0.0038581790 0.0016821844 0.0089322315 -0.0061878037
## 2010-01-08 0.0079225396 -0.0004479241 -0.0066937856 0.0049630301
## 2010-01-11 0.0082095665 -0.0054873060 0.0047825905 0.0132889913
summary = summary(data)
summary
## Index SPY.Adjusted QQQ.Adjusted
## Min. :2010-01-04 00:00:00 Min. : 83.93 Min. : 38.30
## 1st Qu.:2012-03-31 12:00:00 1st Qu.:114.65 1st Qu.: 57.87
## Median :2014-07-02 00:00:00 Median :171.77 Median : 88.29
## Mean :2014-07-01 09:18:41 Mean :167.86 Mean : 91.99
## 3rd Qu.:2016-09-28 12:00:00 3rd Qu.:201.77 3rd Qu.:113.43
## Max. :2018-12-28 00:00:00 Max. :283.49 Max. :183.96
## EEM.Adjusted IWM.Adjusted EFA.Adjusted TLT.Adjusted
## Min. :25.84 Min. : 50.71 Min. :34.29 Min. : 66.00
## 1st Qu.:33.40 1st Qu.: 72.12 1st Qu.:43.81 1st Qu.: 90.83
## Median :35.88 Median :103.45 Median :52.07 Median :102.97
## Mean :36.16 Mean :100.88 Mean :51.33 Mean :100.62
## 3rd Qu.:38.36 3rd Qu.:117.83 3rd Qu.:56.91 3rd Qu.:113.78
## Max. :49.53 Max. :169.05 Max. :70.66 Max. :131.03
## IYR.Adjusted GLD.Adjusted
## Min. :28.89 Min. :100.5
## 1st Qu.:45.38 1st Qu.:116.5
## Median :57.43 Median :122.8
## Mean :56.55 Mean :129.3
## 3rd Qu.:68.75 3rd Qu.:136.5
## Max. :79.10 Max. :184.6
sd = sd(data)
sd
## [1] 49.67862
skewness = apply(data, 2, skewness)
skewness
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## 0.29564618 0.58764146 0.46304528 0.27864230 0.07072858
## TLT.Adjusted IYR.Adjusted GLD.Adjusted
## -0.52173462 -0.16364887 0.98423835
kurtosis = apply(data, 2, kurtosis)
kurtosis
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted EFA.Adjusted
## -1.0076620 -0.6788839 0.2777099 -0.8913111 -0.8925992
## TLT.Adjusted IYR.Adjusted GLD.Adjusted
## -0.7274441 -1.1877635 -0.1329543
#Q2. Based on daily returns and their covariance matrix, and compute weights of minimum variance portfolio (MVP).
Sigma.weekly = cov(dailyreturns)
Sigma.weekly
## SPY QQQ EEM IWM EFA
## SPY 8.824124e-05 9.507893e-05 1.055168e-04 1.049020e-04 9.529041e-05
## QQQ 9.507893e-05 1.196808e-04 1.127970e-04 1.133262e-04 9.891653e-05
## EEM 1.055168e-04 1.127970e-04 1.854857e-04 1.268388e-04 1.357814e-04
## IWM 1.049020e-04 1.133262e-04 1.268388e-04 1.516994e-04 1.143482e-04
## EFA 9.529041e-05 9.891653e-05 1.357814e-04 1.143482e-04 1.331411e-04
## TLT -4.016837e-05 -3.926716e-05 -4.605048e-05 -4.990041e-05 -4.731846e-05
## IYR 7.857693e-05 7.877485e-05 9.944414e-05 1.006593e-04 8.806403e-05
## GLD -1.393566e-07 -1.296188e-06 1.908853e-05 1.967605e-06 9.403897e-06
## TLT IYR GLD
## SPY -4.016837e-05 7.857693e-05 -1.393566e-07
## QQQ -3.926716e-05 7.877485e-05 -1.296188e-06
## EEM -4.605048e-05 9.944414e-05 1.908853e-05
## IWM -4.990041e-05 1.006593e-04 1.967605e-06
## EFA -4.731846e-05 8.806403e-05 9.403897e-06
## TLT 7.704893e-05 -2.202553e-05 1.625844e-05
## IYR -2.202553e-05 1.229320e-04 8.200409e-06
## GLD 1.625844e-05 8.200409e-06 9.977856e-05
ones = rep(1,8)
one.vec = matrix(ones, ncol=1)
a.weekly = inv(Sigma.weekly)%*%one.vec
b.weekly = t(one.vec)%*%a.weekly
mvp.w.weekly =a.weekly / as.numeric(b.weekly)
mvp.w.weekly
## [,1]
## SPY 0.87500522
## QQQ -0.20099862
## EEM -0.12383506
## IWM -0.06301082
## EFA 0.04730848
## TLT 0.43686442
## IYR -0.10171444
## GLD 0.13038081
#Q3. By 2, now use year 2010-2013 data to calculate weekly returns and their covariance matrix, and compute weights of minimum variance portfolio (MVP).
data.weekly <- data[endpoints(data, on="weeks", k=1), ]
#selectByDate(data.weekly,start = "2010/01/01",end = "2013/12/31")
head(data.weekly)
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2010-01-08 93.25535 41.83553 35.05952 55.75293
## 2010-01-15 92.49835 41.20641 34.04507 55.02707
## 2010-01-22 88.89250 39.68758 32.14601 53.34204
## 2010-01-29 87.41111 38.45634 31.06663 51.94217
## 2010-02-05 86.81692 38.62710 30.19015 51.21632
## 2010-02-12 87.94019 39.32809 31.19648 52.72852
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2010-01-08 42.40480 66.67115 31.30751 111.37
## 2010-01-15 42.25665 68.00774 31.11015 110.86
## 2010-01-22 39.90124 68.69465 29.81019 107.17
## 2010-01-29 38.87168 68.92610 29.55837 105.96
## 2010-02-05 38.13099 68.92237 29.65366 104.68
## 2010-02-12 38.33098 67.58110 29.42906 107.04
tail(data.weekly)
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2018-11-23 255.3563 157.1234 38.04298 144.8815
## 2018-11-30 267.3844 167.1502 39.32575 149.5154
## 2018-12-07 255.6667 159.2649 38.17699 141.0903
## 2018-12-14 252.6596 158.9689 38.12913 137.7496
## 2018-12-21 234.8459 145.6359 37.42901 126.1795
## 2018-12-28 241.7244 151.3971 38.12862 130.5929
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2018-11-23 59.16559 111.4575 75.95976 115.77
## 2018-11-30 60.11369 111.6414 78.01066 115.54
## 2018-12-07 58.17917 114.8889 77.98204 118.09
## 2018-12-14 57.90145 114.9568 76.62749 117.06
## 2018-12-21 55.71682 117.4080 71.66856 118.72
## 2018-12-28 56.81064 117.7290 71.68774 121.06
Sigma.weekly = cov(data.weekly)
Sigma.weekly
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## SPY.Adjusted 3051.2655 2165.14274 107.933138 1687.11213
## QQQ.Adjusted 2165.1427 1558.99886 80.643850 1192.86625
## EEM.Adjusted 107.9331 80.64385 16.631783 65.05986
## IWM.Adjusted 1687.1121 1192.86625 65.059859 949.26733
## EFA.Adjusted 431.1944 300.24969 21.861857 245.80679
## TLT.Adjusted 707.8519 493.24495 7.837092 373.82039
## IYR.Adjusted 704.2096 491.19599 19.585912 386.05037
## GLD.Adjusted -540.6809 -365.45347 2.128078 -298.03118
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## SPY.Adjusted 431.19436 707.851905 704.20963 -540.680906
## QQQ.Adjusted 300.24969 493.244952 491.19599 -365.453465
## EEM.Adjusted 21.86186 7.837092 19.58591 2.128078
## IWM.Adjusted 245.80679 373.820389 386.05037 -298.031184
## EFA.Adjusted 70.58884 86.011059 96.77317 -80.968219
## TLT.Adjusted 86.01106 254.229244 190.53520 -75.988040
## IYR.Adjusted 96.77317 190.535195 176.21421 -114.629232
## GLD.Adjusted -80.96822 -75.988040 -114.62923 346.613192
ones = rep(1,8)
one.vec = matrix(ones, ncol=1)
a.weekly = inv(Sigma.weekly)%*%one.vec
b.weekly = t(one.vec)%*%a.weekly
mvp.w.weekly =a.weekly / as.numeric(b.weekly)
mvp.w.weekly
## [,1]
## SPY.Adjusted -0.48556924
## QQQ.Adjusted 0.36413172
## EEM.Adjusted 0.13858224
## IWM.Adjusted 0.05927624
## EFA.Adjusted 0.48189891
## TLT.Adjusted 0.06963643
## IYR.Adjusted 0.42145111
## GLD.Adjusted -0.04940741
#Q4. By 3, now use year 2010-2013 data to calculate monthly returns and their covariance matrix, and compute weights of minimum variance portfolio (MVP).
data.monthly <- data[endpoints(data, on="months", k=1), ]
#selectByDate(data.monthly,start = "2010/01/01",end = "2013/12/31")
head(data.monthly)
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2010-01-29 87.41111 38.45634 31.06663 51.94217
## 2010-02-26 90.13786 40.22682 31.61849 54.26666
## 2010-03-31 95.62543 43.32866 34.18303 58.73316
## 2010-04-30 97.10476 44.30033 34.12624 62.06832
## 2010-05-28 89.38935 41.02549 30.92056 57.39044
## 2010-06-30 84.76425 38.57392 30.48808 52.94647
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2010-01-29 38.87168 68.92610 29.55837 105.96
## 2010-02-26 38.97538 68.69011 31.17139 109.43
## 2010-03-31 41.46411 67.27692 34.21013 108.95
## 2010-04-30 40.30122 69.51173 36.39551 115.36
## 2010-05-28 35.79038 73.06265 34.32696 118.88
## 2010-06-30 35.05241 77.29867 32.72385 121.68
tail(data.monthly)
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## 2018-07-31 271.6643 173.8246 42.94433 162.0657
## 2018-08-31 280.3358 183.8728 41.32650 169.0518
## 2018-09-28 282.0025 183.3550 41.08717 165.1213
## 2018-10-31 262.5150 167.5943 37.48774 146.9780
## 2018-11-30 267.3844 167.1502 39.32575 149.5154
## 2018-12-28 241.7244 151.3971 38.12862 130.5929
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## 2018-07-31 65.96513 114.8179 76.75068 115.99
## 2018-08-31 64.49030 116.3244 78.57381 113.51
## 2018-09-28 65.11279 112.9925 76.33176 112.76
## 2018-10-31 59.81681 109.6812 74.50982 115.15
## 2018-11-30 60.11369 111.6414 78.01066 115.54
## 2018-12-28 56.81064 117.7290 71.68774 121.06
Sigma.monthly = cov(data.monthly)
Sigma.monthly
## SPY.Adjusted QQQ.Adjusted EEM.Adjusted IWM.Adjusted
## SPY.Adjusted 3082.9790 2192.69941 109.773105 1694.4097
## QQQ.Adjusted 2192.6994 1581.88775 82.067073 1200.1213
## EEM.Adjusted 109.7731 82.06707 16.680002 65.5399
## IWM.Adjusted 1694.4097 1200.12131 65.539898 947.9104
## EFA.Adjusted 435.2822 303.75404 22.045643 246.7060
## TLT.Adjusted 706.4301 493.10588 8.017762 371.5877
## IYR.Adjusted 710.3763 497.08651 20.111415 387.5690
## GLD.Adjusted -548.7044 -372.41627 4.103780 -300.7905
## EFA.Adjusted TLT.Adjusted IYR.Adjusted GLD.Adjusted
## SPY.Adjusted 435.28218 706.430090 710.37628 -548.70439
## QQQ.Adjusted 303.75404 493.105881 497.08651 -372.41627
## EEM.Adjusted 22.04564 8.017762 20.11141 4.10378
## IWM.Adjusted 246.70599 371.587732 387.56898 -300.79054
## EFA.Adjusted 71.11492 86.010238 97.51414 -80.28122
## TLT.Adjusted 86.01024 252.456807 189.55029 -79.03892
## IYR.Adjusted 97.51414 189.550292 177.01531 -115.12971
## GLD.Adjusted -80.28122 -79.038918 -115.12971 352.61299
ones = rep(1,8)
one.vec = matrix(ones, ncol=1)
a.monthly = inv(Sigma.monthly)%*%one.vec
b.monthly = t(one.vec)%*%a.monthly
mvp.w.monthly =a.monthly / as.numeric(b.monthly)
mvp.w.monthly
## [,1]
## SPY.Adjusted -0.48141328
## QQQ.Adjusted 0.36139174
## EEM.Adjusted 0.16432491
## IWM.Adjusted 0.05362083
## EFA.Adjusted 0.46451746
## TLT.Adjusted 0.07135831
## IYR.Adjusted 0.42361316
## GLD.Adjusted -0.05741313
#Q5. Compute optimal weights for Q3 and Q4 when short selling is not allowed.
n <- length(tickers)
constraints <- new.constraints(n, lb = -Inf, ub= +Inf)
# sum x.i = 1
constraints <- add.constraints(rep(1,n) , 1 , type = "=", constraints)
ia <- create.historical.ia(data.weekly, 250)
weight <- min.risk.portfolio(ia, constraints)
## Loading required package: kernlab
##
## Attaching package: 'kernlab'
## The following object is masked _by_ '.GlobalEnv':
##
## cross
## The following object is masked from 'package:purrr':
##
## cross
## The following object is masked from 'package:ggplot2':
##
## alpha
weight
## [1] -0.48556924 0.36413172 0.13858224 0.05927624 0.48189891 0.06963643
## [7] 0.42145111 -0.04940741
#Q6. Go to Fama French 3 factors data website:Download Fama/French 3 factor returns’ monthly data (Mkt-RF, SMB and HML).
#a. Compute its mean and covariance matrix using monthly data from 2010-2018.
library(readxl)
FAMA <- read_excel("FAMA.xlsx")
FAMA.ret<-FAMA %>% select(c(2,3,4,5))/100
glimpse(FAMA.ret)
## Observations: 1,125
## Variables: 4
## $ `Mkt-RF` <dbl> 0.0296, 0.0264, 0.0036, -0.0324, 0.0253, 0.0262, -0.000…
## $ SMB <dbl> -0.0230, -0.0140, -0.0132, 0.0004, -0.0020, -0.0004, -0…
## $ HML <dbl> -0.0287, 0.0419, 0.0001, 0.0051, -0.0035, -0.0002, 0.04…
## $ RF <dbl> 0.0022, 0.0025, 0.0023, 0.0032, 0.0031, 0.0028, 0.0025,…
FAMA1 = FAMA %>% select(c(2,3,4,5))
FAMA.mean = apply(FAMA1, 2, mean)
FAMA.cov = cov(FAMA1)
#Q7.Based on CAPM model, compute covariance matrix for the 8-asset portfolio by using past 60 monthly returns from 2014/01 - 2018/12.
Mkt_RF <- FAMA.ret[, 1]
fit3 = lm(formula = cbind(SMB , HML)~Mkt_RF, data=FAMA.ret)
sigF3 = as.matrix(var(cbind(
FAMA.ret$SMB,
FAMA.ret$HML)))
bbeta3 = as.matrix(fit3$coefficients)
bbeta3 = bbeta3[-1,]
bbeta3 = as.matrix(bbeta3)
N <- dim(FAMA.ret)[1]
sigeps3 = crossprod(fit3$residuals)/(N-4)
sigeps3 = diag(diag(sigeps3))
#cov_3f = sigF3*t(bbeta3) *(bbeta3) + sigeps3
#cov_3f
X.3 = cbind(ones, FAMA.ret$Mkt_RF, FAMA.ret$SMB, FAMA.ret$HML)
## Warning in cbind(ones, FAMA.ret$Mkt_RF, FAMA.ret$SMB, FAMA.ret$HML): number
## of rows of result is not a multiple of vector length (arg 1)
b_hat.3 = solve(t(X.3)%*%(X.3))%*%t(X.3)%*%as.matrix(FAMA.ret)
E_hat.3 = as.matrix(FAMA.ret) - X.3%*%b_hat.3
b_hat.3 = as.matrix(b_hat.3[-1,])
diagD_hat.3 = diag(t(E_hat.3)%*%E_hat.3)/(N-4)
#cov_3f.3 = as.matrix(b_hat.3)*sigF3%*%b_hat.3 + diag(diagD_hat.3)
#cov_3f.3