# install.packages("astsa")
# library("astsa")
# install.packages("forecast")
# library(forecast)
# install.packages("tseries")
# library("tseries")
# install.packages("readxl")
library("readxl")

# Uploading and standardizing the dataset
DA = read_excel("DA.xlsx")
y  = as.matrix(DA)[100:nrow(DA),]
dy = apply(log(y),2,diff)
n  = nrow(y)
p  = ncol(y)
m  = apply(y,2,mean)
s  = sqrt(apply(y,2,var))
y  = (y-matrix(m,n,p,byrow=TRUE))%*%diag(1/s)
m  = apply(dy,2,mean)
s  = sqrt(apply(dy,2,var))
dy  = (dy-matrix(m,n-1,p,byrow=TRUE))%*%diag(1/s)

par(mfrow=c(1,1))
ts.plot(y[,1], xlab="Days",ylab="Yield (standardized)",ylim=range(y))
title("Yield Curve - Brazil (jan/06 to dec/16)")
for (i in 2:p)
  lines(y[,i],col=i)
legend("top",legend=names(DA),col=1:8,lty=1)

par(mfrow=c(2,4))
for (i in 1:p)
ts.plot(dy[,i], xlab="Days",ylab="1st diff of yield (standardized)",main=names(DA)[i])


par(mfrow=c(1,1))
ts.plot(dy[,1], xlab="Days",ylab="1st diff of yield (standardized)",ylim=range(dy))
title("Yield Curve - Brazil (jan/06 to dec/16)")
for (i in 2:p)
  lines(dy[,i],col=i)



# Correlations
cor(y)
cor(dy)

# PCA (original data)
pca = princomp(y,scores=TRUE,cor=FALSE)
summary(pca)
prop.var = pca$sdev^2/sum(pca$sdev^2)
pc.y = pca$scores[,1]
pc.y = (pc.y-mean(pc.y))/sqrt(var(pc.y))

par(mfrow=c(1,1))
plot(prop.var,xlab="Number of principal components",ylab="% variance",type="b",ylim=c(0,1))
lines(cumsum(prop.var),type="b",col=2)
abline(h=0,lty=2)
abline(h=1,lty=2,col=2)
legend("right",legend=c("% var","cumulative % var"),col=1:2,lty=1)

par(mfrow=c(1,1))
ts.plot(y[,1], xlab="Days",ylab="Yield (standardized)",ylim=range(y))
title("Yield Curve - Brazil (jan/06 to dec/16)")
for (i in 2:p)
  lines(y[,i],col=i)
legend("top",legend=c(names(DA),"1st PC"),col=c(1:8,1),lty=1,lwd=c(rep(1,8),3))
lines(pc.y,lwd=3)

# PCA (first difference data)
pca = princomp(dy,scores=TRUE,cor=FALSE)
summary(pca)
prop.var = pca$sdev^2/sum(pca$sdev^2)
pc.dy = pca$scores[,1]
pc.dy = (pc.dy-mean(pc.dy))/sqrt(var(pc.dy))

par(mfrow=c(1,1))
plot(prop.var,xlab="Number of principal components",ylab="% variance",type="b",ylim=c(0,1))
lines(cumsum(prop.var),type="b",col=2)
abline(h=0,lty=2)
abline(h=1,lty=2,col=2)
legend("right",legend=c("% var","cumulative % var"),col=1:2,lty=1)

par(mfrow=c(1,1))
ts.plot(dy[,1], xlab="Days",ylab="1st diff of yield (standardized)",ylim=range(dy))
title("Yield Curve - Brazil (jan/06 to dec/16)")
for (i in 2:p)
  lines(dy[,i],col=i)
legend("bottomright",legend=c(names(DA),"1st PC"),col=c(1:8,1),lty=1,lwd=c(rep(1,8),3))
lines(pc.dy,lwd=3)



# Factor analysis
fa.y  = factanal(y,"mle",scores="regression",factors=1)
fa.y
fa.dy = factanal(dy,"mle",scores="regression",factors=1)
fa.dy
cf.y  = fa.y$scores[,1]
cf.dy = fa.dy$scores[,1]
cf.y = (cf.y-mean(cf.y))/sqrt(var(cf.y))
cf.dy = (cf.dy-mean(cf.dy))/sqrt(var(cf.dy))

par(mfrow=c(1,1))
ts.plot(y[,1], xlab="Days",ylab="Yield (standardized)",ylim=range(y))
title("Yield Curve - Brazil (jan/06 to dec/16)")
for (i in 2:p)
  lines(y[,i],col=i)
legend("top",legend=c(names(DA),"1st CF"),col=c(1:8,1),lty=1,lwd=c(rep(1,8),3))
lines(cf.y,lwd=3)

par(mfrow=c(1,1))
ts.plot(dy[,1], xlab="Days",ylab="1st diff of yield (standardized)",ylim=range(dy))
title("Yield Curve - Brazil (jan/06 to dec/16)")
for (i in 2:p)
  lines(dy[,i],col=i)
legend("bottomright",legend=c(names(DA),"1st CF"),col=c(1:8,1),lty=1,lwd=c(rep(1,8),3))
lines(cf.dy,lwd=3)


# Comparing PCA and FA
par(mfrow=c(2,2))
ts.plot(pc.y,xlab="Days",ylab="Index",main="Yield")
lines(cf.y,col=2)
plot(pc.y,cf.y,xlab="Principal component",ylab="Common factor")
abline(0,1,col=2)

ts.plot(pc.dy,xlab="Days",ylab="Index",main="1st diff yield")
lines(cf.dy,col=2)
plot(pc.dy,cf.dy,xlab="Principal component",ylab="Common factor")
abline(0,1,col=2)



# Factor model with stochastic volatility
install.packages("factorstochvol")
library("factorstochvol")
fsv = fsvsample(dy,runningstore=2,runningstoremoments=2)
par(mfrow=c(1,1))
ts.plot(abs(fsv$runningstore$f[1,,1]),xlab="Days",ylab="St.Dev.",main="Common factor")
lines(exp(fsv$runningstore$h[,p+1,1]/2),col=3,lwd=3)

cfsv = fsv$runningstore$f[1,,1]
cfsv = (cfsv-mean(cfsv))/sqrt(var(cfsv))

par(mfrow=c(1,1))
plot(-20+pc.dy,xlab="Days",main="Yield",ylim=c(-30,30),ylab="Common component",axes=FALSE,type="l")
axis(1);box()
lines(cf.dy,col=2)
lines(20+cfsv,col=3)
text(1500,-15,"Principal component",col=1)
text(1500,5,"Commom factor",col=2)
text(1500,25,"Common factor (FSV)",col=3)

pca.load = pca$loadings[,1]/pca$loadings[1,1]
fa.load  = fa.dy$loadings/fa.dy$loadings[1]
fsv.load = t(apply(fsv$facload[,1,],1,quantile,c(0.5,0.05,0.95)))
fsv.load = fsv.load/fsv.load[1,1]
loads = cbind(pca.load,fa.load,fsv.load)
rownames(loads) = names(DA)
colnames(loads) = c("PCA","FA","FSV","5%","95%")
round(loads,3)