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