# Read Ruey S. Tsay's "Analysis of Financial Time Series"

# single-factor model that uses market return as factor
# 13 stock return and sp500 index (market return)
# monthly data from Jan 1990 to Dec 2003


data = read.table("672_2014_factor1.txt", header=T)
n = 168

c = rep(1, n)   # intercept term
X = as.matrix(cbind(c, data[,14]))
Y = as.matrix(data[,1:13])

bet = matrix(0, 2, 13); ehat = matrix(0, n, 13)
for (i in 1:13) {
bet[,i] = solve(t(X)%*%X)%*%(t(X)%*%Y[,i])
ehat[,i] = Y[,i] - X%*%bet[,i]
}

rss = diag(crossprod(ehat))       # RSS
rsq = 1 - rss/diag(var(Y)*(n-1))  # TSS

beta = bet[2,]
result = round(cbind(beta, sqrt(rss/(n-2)), rsq),digits = 3)
colnames(result) <- c("Beta","Sigma", "R-Squared")
result

# industry-factor (BARRA) model that uses industry dummy as beta
# factor is estimated for each period
# 10 stock returns, first 4 financial firms; the next 3 IT firms, and last 3 others
# monthly data from Jan 1990 to Dec 2003


library(Matrix)
data = read.table("672_2014_factor2.txt", header=T)
n = 168
Y = as.matrix(data[,1:10])

d.fin = c(rep(1,4),rep(0,6))
d.it = c(rep(0,4),rep(1,3),rep(0,3))
d.ot = c(rep(0,7),rep(1,3))
X = cbind(d.fin, d.it, d.ot)

f.o = solve(crossprod(X))%*%(t(X)%*%t(Y))  # OLS
ehat.o = t(Y) - X%*%f.o

omega = var(t(ehat.o))
weight = diag(omega^(-1))
w.m = Diagonal(10, x = weight)


x.wls = t(X)%*%w.m%*%X
y.wls = t(X)%*%w.m%*%t(Y)
f.wls = solve(x.wls)%*%y.wls               # WLS

f.1 = ts(f.wls[1,], start=1990, frequency=12)
plot(f.1, main="Factor Realization: financial sector", xlab="Time", type = "l")

f.2 = ts(f.wls[2,], start=1990, frequency=12)
plot(f.2, main="Factor Realization: IT sector", xlab="Time", type = "l")

# statistical-factor model that uses the principle component as factor
# 5 stock returns
# monthly data from Jan 1990 to Dec 1999

data = read.table("672_2014_factor3.txt", header=T)
data = data[1:120,]
mean(data)
v.m = var(data)
e.value = eigen(v.m)$values
e.value
e.vector = eigen(v.m)$vector

pc1 = as.matrix(data)%*%e.vector[,1]
pc2 = as.matrix(data)%*%e.vector[,2]
pc3 = as.matrix(data)%*%e.vector[,3]

var(pc1)
var(pc2)
var(pc3)

cor(cbind(pc1, pc2, pc3))

summary(lm(data[,1]~pc1+pc2+pc3))
summary(lm(data[,5]~pc1+pc2+pc3))