Skip to content

Instantly share code, notes, and snippets.

@tonyday567

tonyday567/rfdata-example.r

Created April 29, 2013 02:25
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save tonyday567/5479381 to your computer and use it in GitHub Desktop.
Save tonyday567/5479381 to your computer and use it in GitHub Desktop.
Download ZIP
example of rates risk factor analysis
Raw
rfdata-example.r
require(plyr)
table <- read.table("/var/folders/h1/yhtz6g614_g999nh_tb4gj600000gn/T/babel-6520FQ/R-import-652nsG",
header=TRUE,
row.names=NULL,
sep="\t",
as.is=TRUE)
specs = table
require(quantmod)
#get US 10y Yield from Fred
getSymbols(specs$code, src="FRED")
usrates = merge.xts(
DGS10 ,
DFII10,
DGS1MO,
DGS1 ,
DGS20 ,
DFII20,
DGS2 ,
DGS30 ,
DFII30,
DGS3MO,
DGS3 ,
DGS5 ,
DFII5 ,
DGS6MO,
DGS7 ,
DFII7)
usrates = usrates / 100
table <- read.table("/var/folders/h1/yhtz6g614_g999nh_tb4gj600000gn/T/babel-6520FQ/R-import-65202M",
header=TRUE,
row.names=NULL,
sep="\t",
as.is=TRUE)
require(quantmod)
getSymbols(table$code, src="FRED")
cpi.monthly = diff(CPIAUCNS)/CPIAUCNS[-length(CPIAUCNS)]
uscash = merge.xts(
DFF/100,
DTB3/100,
DGS1MO/100,
DGS3MO/100
)
colnames(uscash) = c("DFF","DTB3","DGS1MO","DGS3MO")
uscash = uscash["1962::2013"]
# first dates
first.date = function (x) {
index(na.omit(uscash[,x]))[1]
}
start.dates = ldply(1:4,
function (x) {
index(na.omit(uscash[,x]))[1]
})
colnames(start.dates) = "start"
rownames(start.dates) = rownames(uscash)
end.date = index(uscash)[length(uscash)]
ggplot(data=uscash["2001::2013"],aes(x=DTB3,y=DFF)) + geom_point()
ggsave("Fed.TB.png")
uscash$dffxs = uscash$DFF - uscash$DTB3
ggplot(data=uscash["1960::2013"],aes(x=DTB3,y=dffxs)) + geom_point()
uscash$dtb3xs = uscash$DTB3 - uscash$DGS3MO
ggplot(data=uscash["1960::2013"],aes(x=DTB3,y=dtb3xs)) + geom_point()
ggsave("secondary3m.gap.png")
data.fit = na.omit(uscash[,c("dtb3xs","DTB3")])
ss = smooth.spline(x=data.fit$DTB3,y=data.fit$dtb3xs,spar=0.5)
qplot(ss$x,ss$y)
DTB3adj = predict(ss,na.omit(uscash$DTB3))$y
colnames(DTB3adj) = "DTB3adj"
uscash = merge.xts(uscash,DTB3adj=DTB3adj)
# prep dfg (data.frame for graphics)
dfg=data.frame(uscash,Date=as.Date(index(uscash)))
#merging DTB3 with DGS3MO
ggplot(data=dfg,aes(x=Date,y=DTB3))+geom_point()+geom_point(colour="red",aes(x=Date,y=DGS3MO))
ggplot(data=dfg) +
geom_point(colour="red",aes(x=Date,y=DGS3MO)) +
geom_point(colour="blue",aes(x=Date,y=(DTB3 - DTB3adj)))
ggplot(data=dfg) +
geom_point(colour="red",aes(x=Date,y=DGS3MO-(DTB3))) +
geom_point(colour="blue",aes(x=Date,y=DGS3MO-(DTB3-DTB3adj)))
DGS3MOadj = c((uscash$DTB3-uscash$DTB3adj)["1962::1981"],uscash$DGS3MO["1982::2013"])
colnames(DGS3MOadj) = "DGS3MOadj"
uscash = merge.xts(uscash,DGS3MOadj=DGS3MOadj)
dfg=data.frame(uscash,Date=as.Date(index(uscash)))
ggplot(data=dfg) +
geom_point(colour="red",aes(x=Date,y=DGS3MO)) +
geom_point(colour="blue",aes(x=Date,y=DGS3MOadj))
curve.names = c("DGS3MO","DGS6MO","DGS1","DGS2","DGS3","DGS5","DGS7","DGS10","DGS20","DGS30")
duration = aaply(curve.names,1, function (x) {specs[specs$code==x,"duration"]})
order = aaply(curve.names,1, function (x) {which(specs$code==x)})
usnoms = usrates[,order]
xts = merge.xts(usnoms$DGS3MO,DGS3MOadj)
dfg = as.data.frame(xts)
dfg$Date = as.Date(index(xts))
colnames(dfg) = c("DGS3MO","DGS3MOadj","Date")
ggplot(data=dfg) +
geom_point(colour="red",aes(x=Date,y=DGS3MO)) +
geom_point(colour="blue",aes(x=Date,y=DGS3MOadj))
usnoms[,1] = DGS3MOadj[index(usnoms)]
fill.ytm = function (x) {
n = length(x)
last.fill = 0
ytm = rep(NA,n)
for (c in 1:n) {
if (!is.na(x[,c])) {
ytm[c]=x[,c]
if (last.fill == 0) {
ytm[1:(c-1)] = ytm[c]
} else {
ytm[(last.fill+1):(c-1)] = ytm[last.fill] + (duration[(last.fill+1):(c-1)]-duration[last.fill])/(duration[c]-duration[last.fill]) * (ytm[c] - ytm[last.fill])
}
last.fill = c
}
}
fill.ytm=ytm
}
ytm = aaply(usnoms,1,function (x) {fill.ytm(x)})
ytm.xts = xts(ytm,order.by=index(usnoms))
#ytm.xts = ytm.xts[aaply(ytm.xts,1,sum)!=0,]
make.fwd = function (ytm) {
n = length(ytm)
last.fill = 0
fwd = rep(0,n)
for (c in 1:n) {
if (!is.na(ytm[,c])) {
if (last.fill == 0) {
fwd[1:c] = ytm[,c]
} else {
upto = prod((1+fwd[1:last.fill])^(duration[1:last.fill]-c(0,duration[1:(last.fill-1)])))
getto = (1+ytm[,c])^duration[c]
fwd[(last.fill+1):c] = (getto / upto)^(1/(duration[c] - duration[last.fill]))-1
}
last.fill = c
}
}
if (last.fill<n && last.fill>0) {
fwd[(last.fill+1):n]=fwd[last.fill]
}
make.fwd=fwd
}
fwd = aaply(ytm.xts,1,function (x) {make.fwd(x)})
fwd.xts = xts(fwd,order.by=index(usrates))
#fwd.xts = fwd.xts[aaply(fwd.xts,1,sum)!=0,]
fill.ytm.back = function (ytm,fwd) {
n = length(ytm)
last.fill = 0
if (!is.na(ytm[1])) {
for (c in 2:n) {
if (is.na(ytm[c])) {
ytm[c]=prod((1+fwd[1:c])^(duration[1:c]-c(0,duration[1:(c-1)])))^(1/duration[c])-1
}
}
}
fill.ytm.back=ytm
}
ytmf = laply(1:dim(ytm)[1],function (x) {fill.ytm.back(ytm[x,],fwd[x,])})
ytmf.xts = xts(ytmf,order.by=index(usrates))
ytmf.xts = ytmf.xts[aaply(ytmf.xts,1,function (x) {sum(is.na(x))})==0,]
fwd.xts = fwd.xts[aaply(fwd.xts,1,function (x) {sum(x)!=0}),]
png(filename="charts/usytm.png")
require(latticeExtra)
xyplot(ytmf.xts,col=brewer.pal("Blues",n=10)[5:8],
ylim=c(0,0.20),
screens=c(rep(1,4),rep(2,3),rep(3,3)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("3 Month - 2 Year","3-7 Year","10-30 Year"),style=5),
main="US Treasuries since 1962")
dev.off()
png(filename="charts/usytm.recent.png")
require(latticeExtra)
xyplot(ytmf.xts["2007::2013"],col=brewer.pal("Blues",n=10)[5:8],
ylim=c(0,0.05),
screens=c(rep(1,4),rep(2,3),rep(3,3)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("3 Month - 2 Year","3-7 Year","10-30 Year"),style=5),
main="US Treasuries since 2007")
dev.off()
png(filename="charts/usfwd.png")
require(latticeExtra)
xyplot(fwd.xts,col=brewer.pal("Blues",n=10)[5:8],
ylim=c(0,0.2),
screens=c(rep(1,4),rep(2,3),rep(3,3)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("3 Month - 2 Year","3-7 Year","10-30 Year"),style=5),
main="US Forwards since 1962")
dev.off()
png(filename="charts/usytm.recent.png")
require(latticeExtra)
xyplot(fwd.xts["2007::2013"],col=brewer.pal("Blues",n=10)[5:8],
ylim=c(0,0.07),
screens=c(rep(1,4),rep(2,3),rep(3,3)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("3 Month - 2 Year","3-7 Year","10-30 Year"),style=5),
main="US Forwards since 2007")
dev.off()
png(filename="charts/usrp.png")
this = fwd.xts-array(fwd.xts[,3],dim(fwd.xts))
xyplot(this,col=brewer.pal("Blues",n=10)[5:8],
ylim=c(-0.02,0.07),
screens=c(rep(1,4),rep(2,3),rep(3,3)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("3 Month - 2 Year","3-7 Year","10-30 Year"),style=5),
main="US Forwards less 1 year forward")
dev.off()
png(filename="charts/usrp.recent.png")
this = fwd.xts-array(fwd.xts[,3],dim(fwd.xts))
xyplot(this["2007::2013"],col=brewer.pal("Blues",n=10)[5:8],
ylim=c(-0.02,0.07),
screens=c(rep(1,4),rep(2,3),rep(3,3)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("3 Month - 2 Year","3-7 Year","10-30 Year"),style=5),
main="US Forwards less 1 year forward")
dev.off()
curve.names = c("DFII5","DFII7","DFII10","DFII20","DFII30")
duration = aaply(curve.names,1, function (x) {specs[specs$code==x,"duration"]})
order = aaply(curve.names,1, function (x) {which(specs$code==x)})
usreal = usrates[,order]
usreal = usreal[aaply(usreal,1,function (x) {sum(is.na(x))})<dim(usreal)[2],]
start.dates = ldply(1:5,
function (x) {
index(na.omit(usreal[,x]))[1]
})
fill.ytm = function (x) {
n = length(x)
last.fill = 0
ytm = rep(NA,n)
for (c in 1:n) {
if (!is.na(x[,c])) {
ytm[c]=x[,c]
if (last.fill == 0) {
ytm[1:(c-1)] = ytm[c]
} else {
ytm[(last.fill+1):(c-1)] = ytm[last.fill] + (duration[(last.fill+1):(c-1)]-duration[last.fill])/(duration[c]-duration[last.fill]) * (ytm[c] - ytm[last.fill])
}
last.fill = c
}
}
fill.ytm=ytm
}
rytm = aaply(usreal,1,function (x) {fill.ytm(x)})
rytm.xts = xts(rytm,order.by=index(usreal))
#ytm.xts = ytm.xts[aaply(ytm.xts,1,sum)!=0,]
make.fwd = function (ytm) {
n = length(ytm)
last.fill = 0
fwd = rep(0,n)
for (c in 1:n) {
if (!is.na(ytm[,c])) {
if (last.fill == 0) {
fwd[1:c] = ytm[,c]
} else {
upto = prod((1+fwd[1:last.fill])^(duration[1:last.fill]-c(0,duration[1:(last.fill-1)])))
getto = (1+ytm[,c])^duration[c]
fwd[(last.fill+1):c] = (getto / upto)^(1/(duration[c] - duration[last.fill]))-1
}
last.fill = c
}
}
if (last.fill<n && last.fill>0) {
fwd[(last.fill+1):n]=fwd[last.fill]
}
make.fwd=fwd
}
rfwd = aaply(rytm.xts,1,function (x) {make.fwd(x)})
rfwd.xts = xts(rfwd,order.by=index(usreal))
#fwd.xts = fwd.xts[aaply(fwd.xts,1,sum)!=0,]
fill.ytm.back = function (ytm,fwd) {
n = length(ytm)
last.fill = 0
if (!is.na(ytm[1])) {
for (c in 2:n) {
if (is.na(ytm[c])) {
ytm[c]=prod((1+fwd[1:c])^(duration[1:c]-c(0,duration[1:(c-1)])))^(1/duration[c])-1
}
}
}
fill.ytm.back=ytm
}
rytmf = laply(1:dim(rytm)[1],function (x) {fill.ytm.back(rytm[x,],rfwd[x,])})
rytmf.xts = xts(rytmf,order.by=index(usreal))
rytmf.xts = rytmf.xts[aaply(rytmf.xts,1,function (x) {sum(is.na(x))})==0,]
rfwd.xts = rfwd.xts[aaply(rfwd.xts,1,function (x) {sum(x)!=0}),]
png(filename="charts/usrytm.png")
require(latticeExtra)
xyplot(rytmf.xts,col=brewer.pal("Blues",n=10)[5:8],
ylim=c(-0.02,0.05),
screens=c(rep(1,3),rep(2,2)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("5-10 Year","20-30 Year"),style=5),
main="US Real since 2003")
dev.off()
png(filename="charts/usrfwd.png")
require(latticeExtra)
xyplot(rfwd.xts,col=brewer.pal("Blues",n=10)[5:8],
ylim=c(-0.02,0.05),
screens=c(rep(1,3),rep(2,2)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("5-10 Year","20-30 Year"),style=5),
main="US Real since 2003")
dev.off()
breaks = fwd.xts["2003::2013",6:10] - rfwd.xts
str(breaks)
png(filename="charts/usbei.png")
require(latticeExtra)
xyplot(breaks,col=brewer.pal("Blues",n=10)[5:8],
ylim=c(-0.02,0.05),
screens=c(rep(1,3),rep(2,2)),
lattice.options=theEconomist.opts(),
par.settings=theEconomist.theme(box="transparent"),
scale=list(y=list(rot=0)),
strip=strip.custom(factor.levels=c("5-10 Year","20-30 Year"),style=5),
main="US breakeven Inflation since 2003")
dev.off()
cpi.monthly = diff(CPIAUCNS)/CPIAUCNS[-length(CPIAUCNS)]
days.month = function (d) {
library(zoo)
ym <- as.yearmon(d)
days.month=as.numeric(as.Date(ym, frac = 1) - as.Date(ym) + 1)
}
dm = summary(days.month(index(cpi.monthly)))
all.days = index(cpi.monthly)[1]:index(cpi.monthly)[length(cpi.monthly)]
cpi.yearmon = xts(coredata(cpi.monthly)/
days.month(index(cpi.monthly))
,order.by=as.yearmon(index(cpi.monthly)))
cpi.daily = unlist(alply(cpi.yearmon,1,function (x) {
rep(x,days.month(index(x)))
}),use.names=FALSE)
cpi.daily = cpi.yearmon[as.yearmon(as.Date(all.days)),]
cpi = xts(cpi.daily[1:length(all.days)],order.by=as.Date(all.days))
extra = as.Date(index(cpi)[dim(cpi)[1]]+1):as.Date(index(usrates))[dim(usrates)[1]]
extra.xts = xts(rep(cpi[dim(cpi)[1]],length(extra)),order.by=as.Date(extra))
cpi.all = c(cpi,extra.xts)
cpi.daily = cpi.all["2003::2013"]
ryield = lag(array(cpi.daily[index(rfwd.xts)],dim(rfwd.xts))+(rfwd.xts/365))[-1,]*as.integer(diff(index(rfwd.xts)))
by = function (v,m) {
t(t(m)*v)
}
real.return = ryield + by(duration,-diff(rytmf.xts)[-1,])
rrcg5 = 30 * -diff(rytmf.xts)[-1,5]
rrcg = by(duration,-diff(rytmf.xts)[-1,])
which = "2003::2013"
ryield.sum = cumsum(ryield[which])
real.return.sum = cumsum(as.xts(real.return)[which])
rrcg.sum = cumsum(as.xts(rrcg)[which])
dfg = data.frame(ryield=ryield[which],real.return=real.return[which],rrcg=rrcg[which],date=as.Date(index(rytmf.xts[which]))[-1])
dfg.sum = data.frame(rytmf.xts=rytmf.xts[which][-1],rfwd.xts=rfwd.xts[which][-1],ryield=ryield.sum,real.return=real.return.sum,rrcg=rrcg.sum,date=as.Date(index(rytmf.xts[which]))[-1])
dfgm = melt(dfg.sum,id.vars="date",measure.vars=c(
"rytmf.xts.1",
"rytmf.xts.2",
"rytmf.xts.3",
"rytmf.xts.4",
"rytmf.xts.5",
"rfwd.xts.1",
"rfwd.xts.2",
"rfwd.xts.3",
"rfwd.xts.4",
"rfwd.xts.5",
"real.return.1",
"real.return.2",
"real.return.3",
"real.return.4",
"real.return.5"))
dfgm$type = c("fwd","ytm","return")[as.double(is.element(dfgm$variable,as.factor(c(
"real.return.1",
"real.return.2",
"real.return.3",
"real.return.4",
"real.return.5")))) * 2 +
as.double(is.element(dfgm$variable,as.factor(c(
"rytmf.xts.1",
"rytmf.xts.2",
"rytmf.xts.3",
"rytmf.xts.4",
"rytmf.xts.5")))) * 1 +
+1]
ggplot(data=dfgm,aes(x=date,y=value)) +
geom_line(aes(color = variable)) +
facet_grid(type ~ ., scales = "free_y")
p = "2003-01-03::2013-04-25"
nyield = (lag(fwd.xts[,6:10])/365)[-1,]*as.integer(diff(index(fwd.xts)))
iyield = nyield[p] - ryield
nom.cg = by(duration,-diff(ytmf.xts[,6:10])[-1,])
nom.cg.xts = (as.xts(nom.cg))[p]
nom.return = as.xts(as.matrix(nyield[p]) + as.matrix(nom.cg.xts))
ireturn = as.xts(as.matrix(nom.return) - as.matrix(real.return))
iytm = ytmf.xts[p,6:10] - rytmf.xts[p]
ifwd = fwd.xts[p,6:10] - rfwd.xts[p]
ireturn.sum = cumsum(log(1+ireturn))
dfg = data.frame(iytm=iytm,ifwd=ifwd,ireturn.sum=ireturn.sum,date=as.Date(index(iytm)))
dfgm = melt(dfg,id.vars="date",measure.vars=c(
"iytm.6",
"iytm.7",
"iytm.8",
"iytm.9",
"iytm.10",
"ifwd.6",
"ifwd.7",
"ifwd.8",
"ifwd.9",
"ifwd.10",
"ireturn.sum.6",
"ireturn.sum.7",
"ireturn.sum.8",
"ireturn.sum.9",
"ireturn.sum.10"))
dfgm$type = c("fwd","ytm","return")[as.double(is.element(dfgm$variable,as.factor(c(
"ireturn.sum.6",
"ireturn.sum.7",
"ireturn.sum.8",
"ireturn.sum.9",
"ireturn.sum.10")))) * 2 +
as.double(is.element(dfgm$variable,as.factor(c(
"iytm.6",
"iytm.7",
"iytm.8",
"iytm.9",
"iytm.10")))) * 1 +
+1]
ggplot(data=dfgm,aes(x=date,y=value)) +
geom_line(aes(color = variable)) +
facet_grid(type ~ ., scales = "free_y")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment