Macroeconomic Factor Models
options(width = 70, digits=4)
library(ellipse)
##
## Attaching package: 'ellipse'
## The following object is masked from 'package:graphics':
##
## pairs
library(fEcofin) # various data sets
library(PerformanceAnalytics) # performance and risk analysis functions
## 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
library(zoo)
library(readxl)
library(writexl)
data(berndtInvest)
retdata = read_excel("berndt.xlsx")
retdata
## # A tibble: 120 x 18
## date CITCRP CONED CONTIL DATGEN DEC DELTA
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1978-01-31 00:00:00 -0.115 -0.079 -0.129 -0.084 -0.1 -0.028
## 2 1978-02-28 00:00:00 -0.019 -0.003 0.037 -0.097 -0.063 -0.033
## 3 1978-03-31 00:00:00 0.059 0.022 0.003 0.063 0.01 0.07
## 4 1978-04-30 00:00:00 0.127 -0.005 0.18 0.179 0.165 0.15
## 5 1978-05-31 00:00:00 0.005 -0.014 0.061 0.052 0.038 -0.031
## 6 1978-06-30 00:00:00 0.007 0.034 -0.059 -0.023 -0.021 0.023
## 7 1978-07-31 00:00:00 0.032 0.011 0.066 0.143 0.107 0.185
## 8 1978-08-31 00:00:00 0.088 0.024 0.033 0.026 -0.017 -0.021
## 9 1978-09-30 00:00:00 0.011 0.048 -0.013 -0.031 -0.037 -0.081
## 10 1978-10-31 00:00:00 -0.071 -0.067 -0.123 -0.085 -0.077 -0.153
## # ... with 110 more rows, and 11 more variables: GENMIL <dbl>,
## # GERBER <dbl>, IBM <dbl>, MARKET <dbl>, MOBIL <dbl>, PANAM <dbl>,
## # PSNH <dbl>, TANDY <dbl>, TEXACO <dbl>, WEYER <dbl>, RKFREE <dbl>
returns.mat = as.matrix(retdata[, c(-1,-11, -18)])
market.mat = as.matrix(retdata[,11, drop=F])
n.obs = nrow(returns.mat)
X.mat = cbind(rep(1,n.obs),market.mat)
colnames(X.mat)[1] = "intercept"
XX.mat = crossprod(X.mat)
# multivariate least squares
G.hat = solve(XX.mat)%*%crossprod(X.mat,returns.mat)
# can also use solve(qr(X.mat), returns.mat)
beta.hat = G.hat[2,]
E.hat = returns.mat - X.mat%*%G.hat
diagD.hat = diag(crossprod(E.hat)/(n.obs-2))
# compute R2 values from multivariate regression
sumSquares = apply(returns.mat, 2, function(x) {sum( (x - mean(x))^2 )})
R.square = 1 - (n.obs-2)*diagD.hat/sumSquares
R.square
## CITCRP CONED CONTIL DATGEN DEC DELTA GENMIL GERBER
## 0.31777 0.01532 0.11216 0.30363 0.33783 0.12163 0.07919 0.23694
## IBM MOBIL PANAM PSNH TANDY TEXACO WEYER
## 0.27523 0.36882 0.14337 0.01763 0.31986 0.27661 0.43083
# print and plot results
cbind(beta.hat, diagD.hat, R.square)
## beta.hat diagD.hat R.square
## CITCRP 0.66778 0.004511 0.31777
## CONED 0.09102 0.002510 0.01532
## CONTIL 0.73836 0.020334 0.11216
## DATGEN 1.02816 0.011423 0.30363
## DEC 0.84305 0.006564 0.33783
## DELTA 0.48946 0.008152 0.12163
## GENMIL 0.26776 0.003928 0.07919
## GERBER 0.62481 0.005924 0.23694
## IBM 0.45302 0.002546 0.27523
## MOBIL 0.71352 0.004105 0.36882
## PANAM 0.73014 0.015008 0.14337
## PSNH 0.21263 0.011872 0.01763
## TANDY 1.05549 0.011162 0.31986
## TEXACO 0.61328 0.004634 0.27661
## WEYER 0.81687 0.004154 0.43083
par(mfrow=c(1,2))
barplot(beta.hat, horiz=T, main="Beta values", col="blue", cex.names = 0.75, las=1)
barplot(R.square, horiz=T, main="R-square values", col="blue", cex.names = 0.75, las=1)

par(mfrow=c(1,1))
# compute single index model covariance/correlation matrices
cov.si = as.numeric(var(market.mat))*beta.hat%*%t(beta.hat) + diag(diagD.hat)
cor.si = cov2cor(cov.si)
# plot correlations using plotcorr() from ellipse package
rownames(cor.si) = colnames(cor.si)
ord <- order(cor.si[1,])
ordered.cor.si <- cor.si[ord, ord]
plotcorr(ordered.cor.si, col=cm.colors(11)[5*ordered.cor.si + 6])

# compute global min variance portfolio
# use single index covariance
w.gmin.si = solve(cov.si)%*%rep(1,nrow(cov.si))
w.gmin.si = w.gmin.si/sum(w.gmin.si)
colnames(w.gmin.si) = "single.index"
# use sample covariance
w.gmin.sample = solve(var(returns.mat))%*%rep(1,nrow(cov.si))
w.gmin.sample = w.gmin.sample/sum(w.gmin.sample)
colnames(w.gmin.sample) = "sample"
cbind(w.gmin.si, sample = w.gmin.sample)
## single.index sample
## CITCRP 0.043792 -0.060353
## CONED 0.375702 0.376287
## CONTIL 0.005229 -0.002152
## DATGEN -0.023476 -0.065582
## DEC -0.004413 0.036255
## DELTA 0.052499 0.031554
## GENMIL 0.181880 0.197734
## GERBER 0.042721 -0.029664
## IBM 0.186566 0.284569
## MOBIL 0.033722 0.022567
## PANAM 0.007792 0.010707
## PSNH 0.066180 0.075171
## TANDY -0.027191 -0.018677
## TEXACO 0.057823 0.199624
## WEYER 0.001173 -0.058040
par(mfrow=c(2,1))
barplot(t(w.gmin.si), horiz=F, main="Single Index Weights", col="blue", cex.names = 0.75, las=2)
barplot(t(w.gmin.sample), horiz=F, main="Sample Weights", col="blue", cex.names = 0.75, las=2)

par(mfrow=c(1,1))
# compare means and sd values on global min variance portfolios
mu.vals = colMeans(returns.mat)
mu.gmin.si = as.numeric(crossprod(w.gmin.si, mu.vals))
sd.gmin.si = as.numeric(sqrt(t(w.gmin.si)%*%cov.si%*%w.gmin.si))
mu.gmin.sample = as.numeric(crossprod(w.gmin.sample, mu.vals))
sd.gmin.sample = as.numeric(sqrt(t(w.gmin.sample)%*%var(returns.mat)%*%w.gmin.sample))
cbind(mu.gmin.si,mu.gmin.sample, sd.gmin.si, sd.gmin.sample)
## mu.gmin.si mu.gmin.sample sd.gmin.si sd.gmin.sample
## [1,] 0.01365 0.01382 0.03257 0.03264
##
## use lm function to compute single index model regressions for each asset
##
asset.names = colnames(returns.mat)
asset.names
## [1] "CITCRP" "CONED" "CONTIL" "DATGEN" "DEC" "DELTA" "GENMIL"
## [8] "GERBER" "IBM" "MOBIL" "PANAM" "PSNH" "TANDY" "TEXACO"
## [15] "WEYER"
# initialize list object to hold regression objects
reg.list = list()
# loop over all assets and estimate time series regression
i = "CITCRP"
for (i in asset.names) {
reg.df = retdata[, c(i, "MARKET")]
si.formula = as.formula(paste(i,"~", "MARKET", sep=" "))
reg.list[[i]] = lm(si.formula, data=reg.df)
}
# examine the elements of reg.list - they are lm objects!
names(reg.list)
## [1] "CITCRP" "CONED" "CONTIL" "DATGEN" "DEC" "DELTA" "GENMIL"
## [8] "GERBER" "IBM" "MOBIL" "PANAM" "PSNH" "TANDY" "TEXACO"
## [15] "WEYER"
class(reg.list$CITCRP)
## [1] "lm"
reg.list$CITCRP
##
## Call:
## lm(formula = si.formula, data = reg.df)
##
## Coefficients:
## (Intercept) MARKET
## 0.00252 0.66778
summary(reg.list$CITCRP)
##
## Call:
## lm(formula = si.formula, data = reg.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16432 -0.05012 0.00226 0.04351 0.22467
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00252 0.00626 0.40 0.69
## MARKET 0.66778 0.09007 7.41 2e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0672 on 118 degrees of freedom
## Multiple R-squared: 0.318, Adjusted R-squared: 0.312
## F-statistic: 55 on 1 and 118 DF, p-value: 2.03e-11
# plot actual vs. fitted over time
# use chart.TimeSeries() function from PerformanceAnalytics package
dataToPlot = cbind(fitted(reg.list$CITCRP), retdata$CITCRP)
colnames(dataToPlot) = c("Fitted","Actual")
dataToPlot.xts<-as.xts(dataToPlot, order.by = retdata$date)
chart.TimeSeries(dataToPlot.xts, main="Single Index Model for CITCRP",
colorset=c("black","blue"), legend.loc="bottomleft")

# scatterplot of the single index model regression
plot(retdata$MARKET, retdata$CITCRP, main="SI model for CITCRP",
type="p", pch=16, col="blue",
xlab="MARKET", ylab="CITCRP")
abline(h=0, v=0)
abline(reg.list$CITCRP, lwd=2, col="red")

## extract beta values, residual sd's and R2's from list of regression objects
## brute force loop
reg.vals = matrix(0, length(asset.names), 3)
rownames(reg.vals) = asset.names
colnames(reg.vals) = c("beta", "residual.sd", "r.square")
for (i in names(reg.list)) {
tmp.fit = reg.list[[i]]
tmp.summary = summary(tmp.fit)
reg.vals[i, "beta"] = coef(tmp.fit)[2]
reg.vals[i, "residual.sd"] = tmp.summary$sigma
reg.vals[i, "r.square"] = tmp.summary$r.squared
}
reg.vals
## beta residual.sd r.square
## CITCRP 0.66778 0.06716 0.31777
## CONED 0.09102 0.05010 0.01532
## CONTIL 0.73836 0.14260 0.11216
## DATGEN 1.02816 0.10688 0.30363
## DEC 0.84305 0.08102 0.33783
## DELTA 0.48946 0.09029 0.12163
## GENMIL 0.26776 0.06268 0.07919
## GERBER 0.62481 0.07697 0.23694
## IBM 0.45302 0.05046 0.27523
## MOBIL 0.71352 0.06407 0.36882
## PANAM 0.73014 0.12251 0.14337
## PSNH 0.21263 0.10896 0.01763
## TANDY 1.05549 0.10565 0.31986
## TEXACO 0.61328 0.06808 0.27661
## WEYER 0.81687 0.06445 0.43083
# alternatively use R apply function for list objects - lapply or sapply
extractRegVals = function(x) {
# x is an lm object
beta.val = coef(x)[2]
residual.sd.val = summary(x)$sigma
r2.val = summary(x)$r.squared
ret.vals = c(beta.val, residual.sd.val, r2.val)
names(ret.vals) = c("beta", "residual.sd", "r.square")
return(ret.vals)
}
reg.vals = sapply(reg.list, FUN=extractRegVals)
t(reg.vals)
## beta residual.sd r.square
## CITCRP 0.66778 0.06716 0.31777
## CONED 0.09102 0.05010 0.01532
## CONTIL 0.73836 0.14260 0.11216
## DATGEN 1.02816 0.10688 0.30363
## DEC 0.84305 0.08102 0.33783
## DELTA 0.48946 0.09029 0.12163
## GENMIL 0.26776 0.06268 0.07919
## GERBER 0.62481 0.07697 0.23694
## IBM 0.45302 0.05046 0.27523
## MOBIL 0.71352 0.06407 0.36882
## PANAM 0.73014 0.12251 0.14337
## PSNH 0.21263 0.10896 0.01763
## TANDY 1.05549 0.10565 0.31986
## TEXACO 0.61328 0.06808 0.27661
## WEYER 0.81687 0.06445 0.43083
Fundamental Factor Models
# create loading matrix B for industry factor model
n.stocks = ncol(returns.mat)
tech.dum = oil.dum = other.dum = matrix(0,n.stocks,1)
rownames(tech.dum) = rownames(oil.dum) = rownames(other.dum) = asset.names
tech.dum[c(4,5,9,13),] = 1
oil.dum[c(3,6,10,11,14),] = 1
other.dum = 1 - tech.dum - oil.dum
B.mat = cbind(tech.dum,oil.dum,other.dum)
colnames(B.mat) = c("TECH","OIL","OTHER")
# show the factor sensitivity matrix
B.mat
## TECH OIL OTHER
## CITCRP 0 0 1
## CONED 0 0 1
## CONTIL 0 1 0
## DATGEN 1 0 0
## DEC 1 0 0
## DELTA 0 1 0
## GENMIL 0 0 1
## GERBER 0 0 1
## IBM 1 0 0
## MOBIL 0 1 0
## PANAM 0 1 0
## PSNH 0 0 1
## TANDY 1 0 0
## TEXACO 0 1 0
## WEYER 0 0 1
colSums(B.mat)
## TECH OIL OTHER
## 4 5 6
# returns.mat is T x N matrix, and fundamental factor model treats R as N x T.
returns.mat = t(returns.mat)
# Step 1: Estimate OLS F.hat ----
# multivariate OLS regression to estimate K x T matrix of factor returns (K=3)
F.hat = solve(crossprod(B.mat))%*%t(B.mat)%*%returns.mat
# rows of F.hat are time series of estimated industry factors (K X T)
F.hat.df<-data.frame(t(F.hat))
F.hat.df$date<-as.Date(retdata$date)
#
library(reshape2)
library(tidyverse)
## -- Attaching packages ----------------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.1 v purrr 0.3.1
## v tibble 2.0.1 v dplyr 0.8.0.1
## v tidyr 0.8.3 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::first() masks xts::first()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks xts::last()
# plot muliple time series using ggplot
melt(F.hat.df, "date")
## date variable value
## 1 1978-01-31 TECH -0.0720000
## 2 1978-02-28 TECH -0.0517500
## 3 1978-03-31 TECH 0.0335000
## 4 1978-04-30 TECH 0.1322500
## 5 1978-05-31 TECH 0.0620000
## 6 1978-06-30 TECH -0.0155000
## 7 1978-07-31 TECH 0.1340000
## 8 1978-08-31 TECH 0.0700000
## 9 1978-09-30 TECH -0.0547500
## 10 1978-10-31 TECH -0.1035000
## 11 1978-11-30 TECH 0.0562500
## 12 1978-12-31 TECH 0.0860000
## 13 1979-01-31 TECH -0.0080000
## 14 1979-02-28 TECH -0.0582500
## 15 1979-03-31 TECH 0.1082500
## 16 1979-04-30 TECH -0.0210000
## 17 1979-05-31 TECH -0.0607500
## 18 1979-06-30 TECH 0.0347500
## 19 1979-07-31 TECH -0.0232500
## 20 1979-08-31 TECH 0.1210000
## 21 1979-09-30 TECH -0.0237500
## 22 1979-10-31 TECH -0.0930000
## 23 1979-11-30 TECH 0.0967500
## 24 1979-12-31 TECH 0.0157500
## 25 1980-01-31 TECH 0.0782500
## 26 1980-02-29 TECH 0.0395000
## 27 1980-03-31 TECH -0.1257500
## 28 1980-04-30 TECH 0.0050000
## 29 1980-05-31 TECH 0.0635000
## 30 1980-06-30 TECH 0.0577500
## 31 1980-07-31 TECH 0.2380000
## 32 1980-08-31 TECH 0.0807500
## 33 1980-09-30 TECH 0.0112500
## 34 1980-10-31 TECH 0.0007500
## 35 1980-11-30 TECH 0.0757500
## 36 1980-12-31 TECH -0.0112500
## 37 1981-01-31 TECH -0.1277500
## 38 1981-02-28 TECH 0.0067500
## 39 1981-03-31 TECH 0.1322500
## 40 1981-04-30 TECH 0.0487500
## 41 1981-05-31 TECH 0.0930000
## 42 1981-06-30 TECH -0.1180000
## 43 1981-07-31 TECH -0.0012500
## 44 1981-08-31 TECH -0.0802500
## 45 1981-09-30 TECH 0.0110000
## 46 1981-10-31 TECH 0.0972500
## 47 1981-11-30 TECH 0.0165000
## 48 1981-12-31 TECH -0.0362500
## 49 1982-01-31 TECH 0.0507500
## 50 1982-02-28 TECH -0.0790000
## 51 1982-03-31 TECH -0.1227500
## 52 1982-04-30 TECH 0.0572500
## 53 1982-05-31 TECH -0.1070000
## 54 1982-06-30 TECH -0.0310000
## 55 1982-07-31 TECH 0.0057500
## 56 1982-08-31 TECH 0.1482500
## 57 1982-09-30 TECH -0.0527500
## 58 1982-10-31 TECH 0.3010000
## 59 1982-11-30 TECH 0.1362500
## 60 1982-12-31 TECH 0.0180000
## 61 1983-01-31 TECH 0.1377500
## 62 1983-02-28 TECH 0.0672500
## 63 1983-03-31 TECH 0.0245000
## 64 1983-04-30 TECH 0.0632500
## 65 1983-05-31 TECH -0.0542500
## 66 1983-06-30 TECH -0.0135000
## 67 1983-07-31 TECH -0.0227500
## 68 1983-08-31 TECH 0.0032500
## 69 1983-09-30 TECH 0.0160000
## 70 1983-10-31 TECH -0.1230000
## 71 1983-11-30 TECH 0.0167500
## 72 1983-12-31 TECH 0.0590000
## 73 1984-01-31 TECH 0.0175000
## 74 1984-02-29 TECH -0.0307500
## 75 1984-03-31 TECH 0.0347500
## 76 1984-04-30 TECH 0.0650000
## 77 1984-05-31 TECH -0.1100000
## 78 1984-06-30 TECH -0.0107500
## 79 1984-07-31 TECH 0.0015000
## 80 1984-08-31 TECH 0.1692500
## 81 1984-09-30 TECH -0.0632500
## 82 1984-10-31 TECH 0.0140000
## 83 1984-11-30 TECH 0.0072500
## 84 1984-12-31 TECH 0.0455000
## 85 1985-01-31 TECH 0.1442500
## 86 1985-02-28 TECH -0.0417500
## 87 1985-03-31 TECH -0.0560000
## 88 1985-04-30 TECH -0.0860000
## 89 1985-05-31 TECH 0.0237500
## 90 1985-06-30 TECH 0.0020000
## 91 1985-07-31 TECH 0.0165000
## 92 1985-08-31 TECH 0.0152500
## 93 1985-09-30 TECH 0.0007500
## 94 1985-10-31 TECH 0.0545000
## 95 1985-11-30 TECH 0.0720000
## 96 1985-12-31 TECH 0.0822500
## 97 1986-01-31 TECH 0.0102500
## 98 1986-02-28 TECH 0.0370000
## 99 1986-03-31 TECH -0.0372500
## 100 1986-04-30 TECH 0.0382500
## 101 1986-05-31 TECH 0.0245000
## 102 1986-06-30 TECH -0.0872500
## 103 1986-07-31 TECH -0.0450000
## 104 1986-08-31 TECH 0.0862500
## 105 1986-09-30 TECH -0.1222500
## 106 1986-10-31 TECH 0.0922500
## 107 1986-11-30 TECH 0.0497500
## 108 1986-12-31 TECH -0.0372500
## 109 1987-01-31 TECH 0.1807500
## 110 1987-02-28 TECH 0.0917500
## 111 1987-03-31 TECH -0.0192500
## 112 1987-04-30 TECH -0.0042500
## 113 1987-05-31 TECH 0.0147500
## 114 1987-06-30 TECH -0.0107500
## 115 1987-07-31 TECH 0.0337500
## 116 1987-08-31 TECH 0.0572500
## 117 1987-09-30 TECH -0.0002500
## 118 1987-10-31 TECH -0.2640000
## 119 1987-11-30 TECH -0.1197500
## 120 1987-12-31 TECH 0.0995000
## 121 1978-01-31 OIL -0.0464000
## 122 1978-02-28 OIL -0.0192000
## 123 1978-03-31 OIL 0.0642000
## 124 1978-04-30 OIL 0.0992000
## 125 1978-05-31 OIL 0.0144000
## 126 1978-06-30 OIL -0.0170000
## 127 1978-07-31 OIL 0.1050000
## 128 1978-08-31 OIL 0.0198000
## 129 1978-09-30 OIL 0.0110000
## 130 1978-10-31 OIL -0.1322000
## 131 1978-11-30 OIL 0.0218000
## 132 1978-12-31 OIL 0.0092000
## 133 1979-01-31 OIL 0.0078000
## 134 1979-02-28 OIL -0.0444000
## 135 1979-03-31 OIL 0.0562000
## 136 1979-04-30 OIL 0.0332000
## 137 1979-05-31 OIL -0.0236000
## 138 1979-06-30 OIL 0.0514000
## 139 1979-07-31 OIL 0.0820000
## 140 1979-08-31 OIL 0.0332000
## 141 1979-09-30 OIL 0.0234000
## 142 1979-10-31 OIL -0.0912000
## 143 1979-11-30 OIL 0.0328000
## 144 1979-12-31 OIL 0.0402000
## 145 1980-01-31 OIL 0.0362000
## 146 1980-02-29 OIL 0.0476000
## 147 1980-03-31 OIL -0.0992000
## 148 1980-04-30 OIL 0.0514000
## 149 1980-05-31 OIL 0.0596000
## 150 1980-06-30 OIL 0.0110000
## 151 1980-07-31 OIL 0.1168000
## 152 1980-08-31 OIL -0.0290000
## 153 1980-09-30 OIL -0.0280000
## 154 1980-10-31 OIL 0.0404000
## 155 1980-11-30 OIL 0.1642000
## 156 1980-12-31 OIL -0.0394000
## 157 1981-01-31 OIL 0.0086000
## 158 1981-02-28 OIL -0.0116000
## 159 1981-03-31 OIL 0.0172000
## 160 1981-04-30 OIL 0.0068000
## 161 1981-05-31 OIL 0.0394000
## 162 1981-06-30 OIL -0.0310000
## 163 1981-07-31 OIL -0.0672000
## 164 1981-08-31 OIL -0.0450000
## 165 1981-09-30 OIL -0.0640000
## 166 1981-10-31 OIL 0.0020000
## 167 1981-11-30 OIL 0.0496000
## 168 1981-12-31 OIL -0.1038000
## 169 1982-01-31 OIL 0.0376000
## 170 1982-02-28 OIL -0.0184000
## 171 1982-03-31 OIL 0.0082000
## 172 1982-04-30 OIL 0.0206000
## 173 1982-05-31 OIL 0.0154000
## 174 1982-06-30 OIL -0.0296000
## 175 1982-07-31 OIL -0.0954000
## 176 1982-08-31 OIL 0.0604000
## 177 1982-09-30 OIL -0.0540000
## 178 1982-10-31 OIL 0.1730000
## 179 1982-11-30 OIL 0.0458000
## 180 1982-12-31 OIL 0.0456000
## 181 1983-01-31 OIL 0.0566000
## 182 1983-02-28 OIL 0.0830000
## 183 1983-03-31 OIL 0.0296000
## 184 1983-04-30 OIL 0.0486000
## 185 1983-05-31 OIL 0.0102000
## 186 1983-06-30 OIL 0.0628000
## 187 1983-07-31 OIL -0.0504000
## 188 1983-08-31 OIL 0.0530000
## 189 1983-09-30 OIL -0.0276000
## 190 1983-10-31 OIL -0.0120000
## 191 1983-11-30 OIL 0.0626000
## 192 1983-12-31 OIL -0.0058000
## 193 1984-01-31 OIL 0.0352000
## 194 1984-02-29 OIL -0.0444000
## 195 1984-03-31 OIL -0.0272000
## 196 1984-04-30 OIL -0.0900000
## 197 1984-05-31 OIL -0.1932000
## 198 1984-06-30 OIL -0.0144000
## 199 1984-07-31 OIL -0.0640000
## 200 1984-08-31 OIL 0.1054000
## 201 1984-09-30 OIL 0.2080000
## 202 1984-10-31 OIL -0.0374000
## 203 1984-11-30 OIL -0.0220000
## 204 1984-12-31 OIL 0.0398000
## 205 1985-01-31 OIL 0.1294000
## 206 1985-02-28 OIL 0.0112000
## 207 1985-03-31 OIL 0.0254000
## 208 1985-04-30 OIL 0.0078000
## 209 1985-05-31 OIL 0.0804000
## 210 1985-06-30 OIL 0.0052000
## 211 1985-07-31 OIL -0.0004000
## 212 1985-08-31 OIL -0.0044000
## 213 1985-09-30 OIL -0.0288000
## 214 1985-10-31 OIL 0.0580000
## 215 1985-11-30 OIL -0.0198000
## 216 1985-12-31 OIL 0.0188000
## 217 1986-01-31 OIL 0.0384000
## 218 1986-02-28 OIL -0.0134000
## 219 1986-03-31 OIL -0.0136000
## 220 1986-04-30 OIL -0.0304000
## 221 1986-05-31 OIL 0.0244000
## 222 1986-06-30 OIL -0.0710000
## 223 1986-07-31 OIL -0.0282000
## 224 1986-08-31 OIL 0.0616000
## 225 1986-09-30 OIL 0.0122000
## 226 1986-10-31 OIL 0.0440000
## 227 1986-11-30 OIL -0.0244000
## 228 1986-12-31 OIL -0.0320000
## 229 1987-01-31 OIL 0.1498000
## 230 1987-02-28 OIL -0.0370000
## 231 1987-03-31 OIL -0.0184000
## 232 1987-04-30 OIL -0.0266000
## 233 1987-05-31 OIL 0.0524000
## 234 1987-06-30 OIL 0.0724000
## 235 1987-07-31 OIL 0.0576000
## 236 1987-08-31 OIL -0.0440000
## 237 1987-09-30 OIL -0.0498000
## 238 1987-10-31 OIL -0.2548000
## 239 1987-11-30 OIL -0.0562000
## 240 1987-12-31 OIL 0.0390000
## 241 1978-01-31 OTHER -0.0775000
## 242 1978-02-28 OTHER -0.0006667
## 243 1978-03-31 OTHER 0.0220000
## 244 1978-04-30 OTHER 0.0513333
## 245 1978-05-31 OTHER 0.0146667
## 246 1978-06-30 OTHER 0.0190000
## 247 1978-07-31 OTHER 0.0426667
## 248 1978-08-31 OTHER 0.0110000
## 249 1978-09-30 OTHER 0.0318333
## 250 1978-10-31 OTHER -0.0718333
## 251 1978-11-30 OTHER 0.0120000
## 252 1978-12-31 OTHER -0.0238333
## 253 1979-01-31 OTHER 0.0663333
## 254 1979-02-28 OTHER -0.0388333
## 255 1979-03-31 OTHER 0.0180000
## 256 1979-04-30 OTHER -0.0403333
## 257 1979-05-31 OTHER 0.0076667
## 258 1979-06-30 OTHER 0.0583333
## 259 1979-07-31 OTHER -0.0321667
## 260 1979-08-31 OTHER 0.0870000
## 261 1979-09-30 OTHER -0.0188333
## 262 1979-10-31 OTHER -0.0801667
## 263 1979-11-30 OTHER 0.0121667
## 264 1979-12-31 OTHER 0.0205000
## 265 1980-01-31 OTHER 0.0075000
## 266 1980-02-29 OTHER -0.0440000
## 267 1980-03-31 OTHER -0.0781667
## 268 1980-04-30 OTHER 0.0715000
## 269 1980-05-31 OTHER 0.1220000
## 270 1980-06-30 OTHER 0.0180000
## 271 1980-07-31 OTHER 0.0181667
## 272 1980-08-31 OTHER -0.0050000
## 273 1980-09-30 OTHER -0.0238333
## 274 1980-10-31 OTHER -0.0150000
## 275 1980-11-30 OTHER 0.0076667
## 276 1980-12-31 OTHER 0.0765000
## 277 1981-01-31 OTHER 0.0291667
## 278 1981-02-28 OTHER 0.0028333
## 279 1981-03-31 OTHER 0.0768333
## 280 1981-04-30 OTHER 0.0221667
## 281 1981-05-31 OTHER 0.0261667
## 282 1981-06-30 OTHER 0.0385000
## 283 1981-07-31 OTHER -0.0160000
## 284 1981-08-31 OTHER -0.0300000
## 285 1981-09-30 OTHER -0.0151667
## 286 1981-10-31 OTHER 0.0406667
## 287 1981-11-30 OTHER 0.0533333
## 288 1981-12-31 OTHER -0.0441667
## 289 1982-01-31 OTHER -0.0126667
## 290 1982-02-28 OTHER 0.0213333
## 291 1982-03-31 OTHER 0.0321667
## 292 1982-04-30 OTHER 0.0663333
## 293 1982-05-31 OTHER -0.0258333
## 294 1982-06-30 OTHER -0.0121667
## 295 1982-07-31 OTHER -0.0045000
## 296 1982-08-31 OTHER 0.1148333
## 297 1982-09-30 OTHER 0.0516667
## 298 1982-10-31 OTHER 0.1080000
## 299 1982-11-30 OTHER 0.0091667
## 300 1982-12-31 OTHER 0.0003333
## 301 1983-01-31 OTHER 0.0230000
## 302 1983-02-28 OTHER 0.0460000
## 303 1983-03-31 OTHER 0.0571667
## 304 1983-04-30 OTHER 0.0700000
## 305 1983-05-31 OTHER -0.0236667
## 306 1983-06-30 OTHER -0.0046667
## 307 1983-07-31 OTHER -0.0338333
## 308 1983-08-31 OTHER -0.0006667
## 309 1983-09-30 OTHER 0.0235000
## 310 1983-10-31 OTHER 0.0066667
## 311 1983-11-30 OTHER 0.0336667
## 312 1983-12-31 OTHER -0.0300000
## 313 1984-01-31 OTHER 0.0098333
## 314 1984-02-29 OTHER -0.0626667
## 315 1984-03-31 OTHER -0.0255000
## 316 1984-04-30 OTHER -0.0843333
## 317 1984-05-31 OTHER -0.0233333
## 318 1984-06-30 OTHER 0.0378333
## 319 1984-07-31 OTHER -0.0378333
## 320 1984-08-31 OTHER 0.0596667
## 321 1984-09-30 OTHER 0.0390000
## 322 1984-10-31 OTHER -0.0340000
## 323 1984-11-30 OTHER 0.0420000
## 324 1984-12-31 OTHER 0.0051667
## 325 1985-01-31 OTHER 0.1021667
## 326 1985-02-28 OTHER 0.0058333
## 327 1985-03-31 OTHER 0.0215000
## 328 1985-04-30 OTHER -0.0236667
## 329 1985-05-31 OTHER 0.1088333
## 330 1985-06-30 OTHER 0.0560000
## 331 1985-07-31 OTHER 0.0148333
## 332 1985-08-31 OTHER 0.0191667
## 333 1985-09-30 OTHER -0.0340000
## 334 1985-10-31 OTHER 0.0483333
## 335 1985-11-30 OTHER 0.0860000
## 336 1985-12-31 OTHER 0.0493333
## 337 1986-01-31 OTHER 0.0423333
## 338 1986-02-28 OTHER 0.0793333
## 339 1986-03-31 OTHER 0.0823333
## 340 1986-04-30 OTHER -0.0265000
## 341 1986-05-31 OTHER 0.0111667
## 342 1986-06-30 OTHER 0.0090000
## 343 1986-07-31 OTHER -0.0313333
## 344 1986-08-31 OTHER 0.0848333
## 345 1986-09-30 OTHER -0.0896667
## 346 1986-10-31 OTHER 0.0743333
## 347 1986-11-30 OTHER -0.0005000
## 348 1986-12-31 OTHER -0.0318333
## 349 1987-01-31 OTHER 0.1011667
## 350 1987-02-28 OTHER 0.0025000
## 351 1987-03-31 OTHER 0.0013333
## 352 1987-04-30 OTHER -0.0476667
## 353 1987-05-31 OTHER 0.0280000
## 354 1987-06-30 OTHER -0.0123333
## 355 1987-07-31 OTHER 0.0350000
## 356 1987-08-31 OTHER 0.0355000
## 357 1987-09-30 OTHER -0.0786667
## 358 1987-10-31 OTHER -0.1733333
## 359 1987-11-30 OTHER -0.0430000
## 360 1987-12-31 OTHER -0.0026667
F.hat.df %>% melt("date") %>%
ggplot(aes(x = date, y = value, group=variable,color=variable)) +
geom_line() +
scale_x_date()

# Compute residual variance from OLS regression ----
# compute N x T matrix of industry factor model residuals
E.hat = returns.mat - B.mat%*%F.hat
# compute residual variances from time series of errors
diagD.hat = apply(E.hat, 1, var)
Dinv.hat = diag(diagD.hat^(-1))
# Step 2: Run multivariate FGLS regression to estimate K x T matrix of factor returns ----
H.hat = solve(t(B.mat)%*%Dinv.hat%*%B.mat)%*%t(B.mat)%*%Dinv.hat
colnames(H.hat) = asset.names
# note: rows of H sum to one so are weights in factor mimicking portfolios
F.hat.gls = H.hat%*%returns.mat
# show gls factor weights
t(H.hat)
## TECH OIL OTHER
## CITCRP 0.0000 0.00000 0.19918
## CONED 0.0000 0.00000 0.22024
## CONTIL 0.0000 0.09611 0.00000
## DATGEN 0.2197 0.00000 0.00000
## DEC 0.3188 0.00000 0.00000
## DELTA 0.0000 0.22326 0.00000
## GENMIL 0.0000 0.00000 0.22967
## GERBER 0.0000 0.00000 0.12697
## IBM 0.2810 0.00000 0.00000
## MOBIL 0.0000 0.28645 0.00000
## PANAM 0.0000 0.11857 0.00000
## PSNH 0.0000 0.00000 0.06683
## TANDY 0.1806 0.00000 0.00000
## TEXACO 0.0000 0.27561 0.00000
## WEYER 0.0000 0.00000 0.15711
colSums(t(H.hat))
## TECH OIL OTHER
## 1 1 1
# compare OLS and GLS fits
F.hat.gls.zoo = zoo(t(F.hat.gls), as.Date(retdata$date))
par(mfrow=c(3,1))
plot(merge(F.hat.gls.zoo[,1], F.hat.gls.zoo[,1]), plot.type="single",
main = "OLS and GLS estimates of TECH factor",
col=c("black", "blue"), lwd=2, ylab="Return")
legend(x = "bottomleft", legend=c("OLS", "GLS"), col=c("black", "blue"), lwd=2)
abline(h=0)
plot(merge(F.hat.gls.zoo[,2], F.hat.gls.zoo[,2]), plot.type="single",
main = "OLS and GLS estimates of OIL factor",
col=c("black", "blue"), lwd=2, ylab="Return")
legend(x = "bottomleft", legend=c("OLS", "GLS"), col=c("black", "blue"), lwd=2)
abline(h=0)
plot(merge(F.hat.gls.zoo[,3], F.hat.gls.zoo[,3]), plot.type="single",
main = "OLS and GLS estimates of OTHER factor",
col=c("black", "blue"), lwd=2, ylab="Return")
legend(x = "bottomleft", legend=c("OLS", "GLS"), col=c("black", "blue"), lwd=2)
abline(h=0)

par(mfrow=c(1,1))
# compute sample covariance matrix of estimated factors
cov.ind = B.mat%*%var(t(F.hat.gls))%*%t(B.mat) + diag(diagD.hat)
cor.ind = cov2cor(cov.ind)
# plot correlations using plotcorr() from ellipse package
rownames(cor.ind) = colnames(cor.ind)
ord <- order(cor.ind[1,])
ordered.cor.ind <- cor.ind[ord, ord]
plotcorr(ordered.cor.ind, col=cm.colors(11)[5*ordered.cor.ind + 6])

# compute industry factor model R-square values
r.square.ind = 1 - diagD.hat/diag(cov.ind)
ind.fm.vals = cbind(B.mat, sqrt(diag(cov.ind)), sqrt(diagD.hat), r.square.ind)
colnames(ind.fm.vals) = c(colnames(B.mat), "fm.sd", "residual.sd", "r.square")
ind.fm.vals
## TECH OIL OTHER fm.sd residual.sd r.square
## CITCRP 0 0 1 0.07291 0.05468 0.4375
## CONED 0 0 1 0.07092 0.05200 0.4624
## CONTIL 0 1 0 0.13258 0.11807 0.2069
## DATGEN 1 0 0 0.10646 0.07189 0.5439
## DEC 1 0 0 0.09862 0.05968 0.6338
## DELTA 0 1 0 0.09817 0.07747 0.3773
## GENMIL 0 0 1 0.07013 0.05092 0.4728
## GERBER 0 0 1 0.08376 0.06849 0.3315
## IBM 1 0 0 0.10102 0.06356 0.6041
## MOBIL 0 1 0 0.09118 0.06839 0.4374
## PANAM 0 1 0 0.12222 0.10630 0.2435
## PSNH 0 0 1 0.10601 0.09440 0.2069
## TANDY 1 0 0 0.11159 0.07930 0.4950
## TEXACO 0 1 0 0.09218 0.06972 0.4279
## WEYER 0 0 1 0.07821 0.06157 0.3802
# compute global minimum variance portfolio
w.gmin.ind = solve(cov.ind)%*%rep(1,nrow(cov.ind))
w.gmin.ind = w.gmin.ind/sum(w.gmin.ind)
t(w.gmin.ind)
## CITCRP CONED CONTIL DATGEN DEC DELTA GENMIL GERBER
## [1,] 0.1401 0.1549 0.02605 0.005638 0.008182 0.0605 0.1615 0.0893
## IBM MOBIL PANAM PSNH TANDY TEXACO WEYER
## [1,] 0.007214 0.07763 0.03213 0.047 0.004635 0.07469 0.1105
# compare weights with weights from sample covariance matrix
par(mfrow=c(2,1))
barplot(t(w.gmin.ind), horiz=F, main="Industry FM Weights", col="blue", cex.names = 0.75, las=2)
barplot(t(w.gmin.sample), horiz=F, main="Sample Weights", col="blue", cex.names = 0.75, las=2)

par(mfrow=c(1,1))
# compare means and sd values on global min variance portfolios
mu.gmin.ind = as.numeric(crossprod(w.gmin.ind, mu.vals))
sd.gmin.ind = as.numeric(sqrt(t(w.gmin.ind)%*%cov.ind%*%w.gmin.ind))
cbind(mu.gmin.sample,mu.gmin.sample, sd.gmin.ind, sd.gmin.sample)
## mu.gmin.sample mu.gmin.sample sd.gmin.ind sd.gmin.sample
## [1,] 0.01382 0.01382 0.05043 0.03264