# FactorExamples.ssc			script file for examples used in Factor model chapter
#
# author: Eric Zivot
# created: January 25, 2002
# revised: May 16, 2002
#

#
# single index model estimation using Berndt data
#
# return data are in matrix returns
# market return data is in vector market

colIds(berndt.dat)
returns = as.matrix(seriesData(berndt.dat[, c(-10, -17)]))
market = as.vector(seriesData(berndt.dat)[,10])
n.obs = nrow(returns)
X.mat = cbind(rep(1,n.obs),market)
G.hat = solve(X.mat,returns)
beta.hat = G.hat[2,]
E.hat = returns - X.mat%*%G.hat
diagD.hat = diag(crossprod(E.hat)/(n.obs-2))
names(diagD.hat) = colIds(G.hat)
r.square = 1 - (n.obs-2)*diagD.hat/diag(var(returns,SumSquares=T))

# print results
t(rbind(beta.hat,sqrt(diagD.hat),r.square))
# plot results
par(mfrow=c(1,2))
barplot(beta.hat,names=names(beta.hat),horiz=T,
main="Beta values")
barplot(r.square,names=names(r.square),horiz=T,
main="R-square values")

#
# compute covariance/correlation matrices
#
cov.si = var(market)*(beta.hat%o%beta.hat) + diag(diagD.hat)
sd = sqrt(diag(cov.si))
cor.si = cov.si/outer(sd,sd)
print(cor.si,digits=1,width=2)
print(cor(returns),digits=1,width=2)

# compute global min variance portfolio
w.gmin.si = solve(cov.si)%*%rep(1,nrow(cov.si))
w.gmin.si = w.gmin.si/sum(w.gmin.si)
t(w.gmin.si)

w.gmin.sample = solve(var(returns))%*%rep(1,nrow(cov.si))
w.gmin.sample = w.gmin.sample/sum(w.gmin.sample)
t(w.gmin.sample)

#
# macroeconomic factor model using Berndt data
#
infl = getReturns(CPI.dat)
ipg = getReturns(IP.dat)

# estimate VAR(1) over period Jan 1978 - Dec 1987)
factor.ts = seriesMerge(ipg,infl)
var6.fit = VAR(cbind(CPI,IP)~ar(6),data=factor.ts,
start="July 1977",end="Jan 1988",in.format="%m %Y")
factor.surprise = residuals(var6.fit)

returns = as.matrix(seriesData(berndt.dat[, c(-10, -17)]))
factor.surprise = as.matrix(seriesData(factor.surprise))
n.obs = nrow(returns)
X.mat = cbind(rep(1,n.obs),factor.surprise)
G.hat = solve(X.mat,returns)
beta.hat = t(G.hat[2:3,])
E.hat = returns - X.mat%*%G.hat
diagD.hat = diag(crossprod(E.hat)/(n.obs-3))
names(diagD.hat) = colIds(G.hat)
r.square = 1 - (n.obs-2)*diagD.hat/diag(var(returns,SumSquares=T))

par(mfrow=c(1,3))
barplot(beta.hat[,1],names=names(beta.hat),horiz=T,
main="Beta values for inflation surprise")
barplot(beta.hat[,2],names=names(beta.hat),horiz=T,
main="Beta values for IP growth surprise")
barplot(r.square,names=names(r.square),horiz=T,
main="R-square values")
par(mfrow=c(1,1))

# compute covariance
cov.macro = beta.hat%*%var(factor.surprise)%*%t(beta.hat) +
diag(diagD.hat)
sd = sqrt(diag(cov.macro))
cor.macro = cov.macro/outer(sd,sd)
print(cor.macro,digits=1,width=2)

w.gmin.macro = solve(cov.macro)%*%rep(1,nrow(cov.macro))
w.gmin.macro = w.gmin.macro/sum(w.gmin.macro)
t(w.gmin.macro)

#
# Estimate industry factor model for Berndt data
#

returns = as.matrix(seriesData(berndt.dat[, c(-10, -17)]))
# create loading matrix B for industry factor model

n.stocks = numCols(returns)
tech.dum = oil.dum = other.dum = matrix(0,n.stocks,1)
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 = cbind(tech.dum,oil.dum,other.dum)
dimnames(B) = list(colIds(returns),c("TECH","OIL","OTHER"))
B

# returns are in T x N matrix. 
# Note that Factor model treats R as N x T.
returns = t(returns)
# multivariate OLS regression to estimate K x T matrix of factor returns 
F.hat = solve(crossprod(B))%*%t(B)%*%returns
# compute N x T matrix of industry factor model residuals
E.hat = returns - B%*%F.hat
# compute residual variances from time series of errors
diagD.hat = rowVars(E.hat)
Dinv.hat = diag(diagD.hat^(-1))
# multivariate FGLS regression to estimate K x T matrix of factor returns 
H = solve(t(B)%*%Dinv.hat%*%B)%*%t(B)%*%Dinv.hat
t(H)
# note: rows of H sum to one
rowSums(H)
# create factor mimicking portfolios
F.hat = H%*%returns
F.hat = t(F.hat)	
tsplot(F.hat)
legend(0,-0.15,legend=colIds(F.hat),lty=1:3)
# compute weights in factor mimicking portfolio: not needed since
# rows of H already sum to 1
# weights = t(H/rowSums(H))
# row.names(weights) = colIds(returns)

# compute sample covariance matrix of estimated factors

cov.ind = B%*%var(F.hat)%*%t(B) + diag(diagD.hat)
sd = sqrt(diag(cov.ind))
cor.ind = cov.ind/outer(sd,sd)
print(cor.ind,digits=1,width=2)

# 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)

#
# compute Factor model decomposition of Berndt data
#

returns = as.matrix(seriesData(berndt.dat[, c(-10, -17)]))
args(factanal)
factor.fit = factanal(returns,factors=2,method="mle")
factor.fit = factanal(returns,factors=3,method="mle")

# method functions: biplot, fitted, loadings, print, residuals, 
# rotate, summary
factor.fit
summary(factor.fit)
plot(loadings(factor.fit))

# rotate factors
factor.fit2 = rotate(factor.fit,rotation="quartimax")
loadings(factor.fit2)

fitted(factor.fit)		# factor analysis correlation matrix
factor.ret = predict(factor.fit,type="weighted.ls")


# extract properly scaled loadings and residual variances
S.hat = diag(factor.fit$scale)
D.hat = S.hat%*%diag(factor.fit$uniqueness)%*%S.hat
D.hat.inv = diag(1/diag(D.hat))
B.hat = S.hat%*%loadings(factor.fit)
# compute factor analysis covariance
cov.fa = B.hat%*%t(B.hat)+D.hat
dimnames(cov.fa) = list(colIds(returns),colIds(returns))
# compute factor mimicking portfolios
W.hat = solve(t(B.hat)%*%D.hat.inv%*%B.hat)%*%t(B.hat)%*%D.hat.inv
W.hat = W.hat/rowSums(W.hat)
FMP.hat = returns%*%t(W.hat)

# global minimum variance portfolio

w.gmin.fa = solve(cov.fa)%*%rep(1,nrow(cov.fa))
w.gmin.fa = w.gmin.fa/sum(w.gmin.fa)
t(w.gmin.fa)


par(mfrow=c(1,2))
barplot(W.hat[1,],names=colIds(returns),horiz=T,
main="Beta values for first factor")
barplot(W.hat[2,],names=colIds(returns),horiz=T,
main="Beta values for second factor")
par(mfrow=c(1,1))

#
# principal component analysis of Berndt dta
#
args(mfactor)
returns.ts = berndt.dat[,c(-10,-17)]
returns.ts = timeSeries(pos=positions(returns.ts),
data=as.matrix(seriesData(returns.ts)))

pc.mfactor = mfactor(returns.ts)
class(pc.mfactor)
names(pc.mfactor)
# see help for mfactor.object
# methods:
# factors, loadings, plot, print, residuals, vcov, mfactor.r2
loadings(pc.mfactor)
mfactor.r2(pc.mfactor)
plot(pc.mfactor, ask=F)

cov.pca = vcov(pc.mfactor)
sd = sqrt(diag(cov.pca))
cor.pca = cov.pca/outer(sd,sd)
print(cor.pca,digits=1,width=2)

# global min variance portfolio
w.gmin.pca = solve(cov.pca)%*%rep(1,nrow(cov.pca))
w.gmin.pca = w.gmin.pca/sum(w.gmin.pca)
t(w.gmin.pca)


pc2.mfactor = mfactor(returns.ts,k=2)
pc2.mfactor
plot(pc2.mfactor, ask=F)
# plot of factor returns
fplot(factors(pc2.mfactor))

# specialized functions:
# mfactor.r2, mimic
# note: mimic doesn't work if the data in the data slot
# is a data.frame.

pc2.mimic = mimic(pc2.mfactor)
class(pc2.mimic)
pc2.mimic.sum = summary(pc2.mimic,n.top=3)
pc2.mimic.sum
par(mfrow=c(1,2))
plot(pc2.mimic.sum)

# plot loadings on factors

pc2.betas = loadings(pc2.mfactor)
par(mfrow=c(1,2))
barplot(pc2.betas[1,],names=colIds(pc2.betas),horiz=T,
main="Beta values for first PCA factor")
barplot(pc2.betas[2,],names=colIds(pc2.betas),horiz=T,
main="Beta values for second PCA factor")
par(mfrow=c(1,1))

#
# APA
#
folio.mf = mfactor(folio.dat,k=15)
folio.mf
screeplot.mfactor(folio.mf)
fplot(factors(folio.mf))

folio.m = mimic(folio.mf)
folio.ms = summary(folio.m,n.top=5)
folio.ms

folio.cov = vcov(folio.mf)
sd = sqrt(diag(folio.cov))
folio.cor = folio.cov/outer(sd,sd)
top.names = c(as.character(folio.ms[[1]][,1]),
rev(as.character(folio.ms[[1]][,3])))
image.plot(folio.cor[top.names, top.names],
sub="Risk factor 1", main="Correlations of top positions")

# determine factors using CK method
folio.mf.ck = mfactor(folio.dat,k="ck",max.k=10,sig=0.05)
folio.mf.bn = mfactor(folio.dat,k="bn",max.k=10,sig=0.05)