######################################################################
###########################Pricing Functions##########################
######################################################################

#Issues:
#payment.dates

#Conventions:
#Today is 1 day before cash maturity
#Ex-coupon date at 7 => need to change back

#p(...) gives price as at date of settlement

#a+b in mod m
mod<-function(a,b,m){
res<-round( ((a+b)/m - floor((a+b)/m))*m )
ifelse(res==0,m,res)}

#next interest payment date: date = date, maturity = maturity of bond
nipd<-function(date,maturity){
y<-as.numeric(format.Date(date,"%Y"))
m.m<-as.numeric(format.Date(maturity,"%m"))
d.m<-as.numeric(format.Date(maturity,"%d"))
p<-as.numeric(c(as.Date(ISOdate(y,m.m,d.m)),as.Date(ISOdate(y,mod(m.m,6,12),d.m)),as.Date(ISOdate(y+1,m.m,d.m)),as.Date(ISOdate(y+1,mod(m.m,6,12),d.m)))-date)
for(i in 1:4){if(p[i]<0 || is.na(p[i]))(p[i]<-100000000)}	#|| is.na(p[i]) added to work with R 2.3
date+min(p)}

#days in half year: date is output from nipd
dihy<-function(date){
y<-as.numeric(format.Date(date,"%Y"))
m<-as.numeric(format.Date(date,"%m"))
d<-as.numeric(format.Date(date,"%d"))
ifelse(m>6,as.numeric(date-as.Date(ISOdate(y,m-6,d))),as.numeric(date-as.Date(ISOdate(y-1,m+6,d))))}

#bond price: y is ytm, c is coupon rate, s is days till settlement
p1<-function(date,maturity,y,c,s.b){
i<-y/200
g<-c/2
v<-1/(1+i)
#dos = date of settlement
#nip = next interest payment
dos<-date+s.b
nip<-nipd(date,maturity)
f<-as.numeric(nip-dos)
d<-dihy(nip)
n<-round(as.numeric(maturity-nipd(date,maturity))/182.5,0)
a<-(1-v^n)/i
v^(f/d)*(g*(1+a)+100*v^n)}

#bond price ex interest: y is ytm, c is coupon rate, s.b is days till settlement
p2<-function(date,maturity,y,c,s.b){
i<-y/200
g<-c/2
v<-1/(1+i)
#dos = date of settlement
#nip = next interest payment
dos<-date+s.b
nip<-nipd(date,maturity)
f<-as.numeric(nip-dos)
d<-dihy(nip)
n<-round(as.numeric(maturity-nipd(date,maturity))/182.5,0)
a<-(1-v^n)/i
v^(f/d)*(g*a+100*v^n)}

#OIS equivalent price
p3<-function(date,maturity,y,s.o){
f<-as.numeric(maturity-(date+s.o))
100/(1+(f/365)*y/100)}

#general price funtion: y is ytm, c is coupon rate, s is days till settlement, s.o is days till OIS settlement
#1. c=0 then p3, else
#2. as.numeric(nipd(date,maturity)-dos)<(8+s) then p2, else p1
p<-function(date,maturity,y,c,s.b,s.o){
ifelse(c==0,p3(date,maturity,y,s.o),ifelse(as.numeric(nipd(date,maturity)-(date+s.b))<8,p2(date,maturity,y,c,s.b),p1(date,maturity,y,c,s.b)))}

#price to yield
pty<-function(date,maturity,c,price,s.b,s.o,rf){
y<-0
y[1]<-rf-6
y[3]<-rf+6
y[2]<-rf
for(i in 1:15){
pr<-p(date,maturity,y[2],c,s.b,s.o)
ifelse(pr>=price,{y[1]<-y[2];y[2]<-(y[2]+y[3])/2},{y[3]<-y[2];y[2]<-(y[1]+y[2])/2})}
round(y[2],3)}

#pty for whole data set
make.y<-function(dates,Coupon,Phat,s.b,s.o,rf){
y<-0
n<-length(dates)
for(i in 2:n){
y[i-1]<-pty(dates[1]-1,dates[i],Coupon[i],Phat[i-1],s.b,s.o,rf)}
y}

#Remove bonds with maturity less than 1.5 years.
clean<-function(dat){
b1<-dat[,1]
b2<-dat[,2]*100
b3<-as.Date(dat[,3])
b4<-dat[,4]
#the list of offendors
	#First where do we cut them off: before 13/10/2000 at 1 year, after at 1.5 years
	KillNumber<-584
	ifelse(b3[1]<=as.Date(ISOdate(2000,10,13)),KillNumber<-365,KillNumber<-548)
list<-!{{(b3-b3[1]+1)<KillNumber}&{b2!=0}}
#Cleaning
j<-0
Issue.ID<-b1[1]
Coupon<-b2[1]
Maturity<-b3[1]
Yield<-b4[1]
for(i in 1:length(b1)){
j<-(j+1)
ifelse(list[i],{Issue.ID[j]<-b1[i];Coupon[j]<-b2[i];Maturity[j]<-b3[i];Yield[j]<-b4[i]},j<-(j-1))}
data.frame(Issue.ID,Coupon,Maturity,Yield)}

#Which settlement dates (ignoring public holidays)
set<-function(date){
s<-2
s[2]<-1
if(weekdays(date,T)=="Thu"){s[1]<-4}
if(weekdays(date,T)=="Fri"){s[1]<-4;s[2]<-3}
if(weekdays(date,T)=="Sat"){s[1]<-3;s[2]<-2}
s}

#imputs a date and outputs last week day at or before that date
tradeday<-function(x){
ifelse(weekdays(x)=="Sunday",x<-(x-2),ifelse(weekdays(x)=="Saturday",x<-(x-1),x))
x}

######################################################################
###########################Discount Function##########################
######################################################################

#Date function: So 2005-02-31 becomes 2005-02-28 not 2005-03-03
ISOdate2<-function(y,m,d){
m.o<-20
while(abs(m-m.o)>0.1){
o<-ISOdate(y,m,d)
m.o<-as.numeric(format.Date(o,"%m"))
d<-d-1}
o}

#The dates of payment of coupons - a sub function: from today
payment.dates<-function(date,maturity){
dates<-c(date,nipd(date,maturity))
d<-as.numeric(format.Date(dates[2],"%d"))
i<-2
while(abs(as.numeric(max(dates)-maturity))>6){
	y<-as.numeric(format.Date(dates[i],"%Y"))
	m<-as.numeric(format.Date(dates[i],"%m"))
	dates<-c(dates,as.Date(ISOdate2(ifelse(m<7,y,y+1),mod(m,6,12),d)))
	i<-i+1}
dates}

#The times to payment of coupons from TODAY
tau<-function(date,maturity,c,s.b){
pd<-payment.dates(date,maturity)
tau.1<-0
for(i in 2:length(pd)){tau.1[i-1]<-as.numeric(pd[i]-pd[1])}
tau.2<-0
if(length(tau.1)==1)(tau.2<-tau.1)
if(length(tau.1)==2)(tau.2<-tau.1[2])
if(length(tau.1)>2)(tau.2<-tau.1[2:length(tau.1)])
x<-0
ifelse(c==0, x<-as.numeric(maturity-date),
ifelse(tau.1[1]<(8+s.b),x<-tau.2,x<-tau.1))
x
}

#coupon stream
coupon<-function(date,maturity,c,s.b){
l<-length(tau(date,maturity,c,s.b))
ifelse(l==1,x<-100+c/2,x<-c(rep(c/2,l-1),100+c/2))
x}

#Duration of a security
dur<-function(date,maturity,c,s.b){
co<-coupon(date,maturity,c,s.b)
ct<-sum(co)
t<-tau(date,maturity,c,s.b)/365
sum(t*co/ct)}

######################################################################
###########################MLES Functions#############################
######################################################################

#f_k(t)
f<-function(t,k,a){1/(1+k*a*t)}

#Make H
make.H<-function(dat,a,s.b,s.o,D){
N<-dim(dat)[1]
H<-matrix(0,N-1,D)
attach(dat, warn.conflicts=FALSE)
dates<-as.Date(dat[,3])
date<-(dates[1]-1)
for(i in 2:N){
	c<-Coupon[i]
	ifelse(c==0,s<-s.o,s<-s.b)	#settlement date
	maturity<-dates[i]
	c_ij<-coupon(date,maturity,c,s)
	t_ij<-tau(date,maturity,c,s)/365
	for(k in 1:D){
		f_ij<-(f(t_ij,k,a)/f(s/365,k,a))	#rescale basis functions
		H[i-1,k]<-sum(c_ij*f_ij)
}}
H}

#make P vector - NOT including cash
make.P<-function(dat,s.b,s.o){
attach(dat, warn.conflicts=FALSE)
dates<-as.Date(dat[,3])
N<-dim(dat)[1]-1
P<-0
for(i in 1:N){P[i]<-p(dates[1]-1,dates[i+1],Yield[i+1],Coupon[i+1],s.b,s.o)}
P}

#make W vector: 1/duration - not including cash
make.W<-function(dat,s.b){
dates<-as.Date(dat[,3])
tmp<-0
for(i in 2:length(dates)){tmp[i-1]<-dur(dates[1]-1,dates[i],dat[i,2],s.b)}
diag(1/tmp)}

#plot the discount function
d<-function(Zhat,l,a,jdi){
idx<-seq(0,l,jdi)
dfv<-0
for(i in 1:length(idx)){
tmp<-0
for(j in 1:length(Zhat)){tmp<-tmp+Zhat[j]*f(idx[i],j,a)}
dfv[i]<-tmp}
data.frame(idx,dfv)}

#plot the zero coupon yield
yc<-function(dc,y){
for(i in 2:dim(dc)[1]){y[i]<-100*(-log(dc[i,2])/dc[i,1])}
data.frame(dc[,1],y)}

#discount function for use in inst forward, offset by 0.005 years
d2<-function(Zhat,l,a,jdi){
idx<-seq(0,l,jdi)+0.005
dfv<-0
for(i in 1:length(idx)){
tmp<-0
for(j in 1:length(Zhat)){tmp<-tmp+Zhat[j]*f(idx[i],j,a)}
dfv[i]<-tmp}
data.frame(idx,dfv)}

# Forward
fc<-function(yc,yc2,rf){
y.c<-yc[,2]/100
y2.c<-yc2[,2]/100
t<-yc[,1]
t2<-yc2[,1]
n<-length(t)
f.c<-c(rf,(1/(t2[2]-t[2]))*100*(t2[2:n]*y2.c[2:n] - t[2:n]*y.c[2:n]))
data.frame(t,f.c)}

#Used to generate (and then plot) the basis functions
basis<-function(a,l,j,d){
idx<-seq(0,l,j)
n<-length(idx)
tmp<-matrix(0,n,(d+1))
tmp[,1]<-idx
for(i in 1:d){tmp[,(i+1)]<-f(idx,i,a)}
tmp}

######################################################################
###########################go<-function###############################
######################################################################

#s.b	days to settlement of bonds
#s.o	days to settlement of OIS
#a.i	scalling factor
#d.i	number of basis ellements
#j.i	interval length of plotting graphs
#l.i	number of years to plot graph out to

#s.b=-1;s.o=-1;a.i=0.05;d.i=8;j.i=0.02;l.i=10;graph=FALSE

go<-function(dat,s.b=-1,s.o=-1,a.i=0.05,d.i=8,j.i=0.02,l.i=10){
dat<-clean(dat)
if(s.b==-1){s<-set(as.Date(dat[1,3])-1);s.b<-s[1]}
if(s.o==-1){s<-set(as.Date(dat[1,3])-1);s.o<-s[2]}
attach(dat, warn.conflicts=FALSE)
dates<-Maturity
#risk free rate
rf<-Yield[1]
#major matricies
H<-make.H(dat,a.i,s.b,s.o,d.i)
P<-make.P(dat,s.b,s.o)
W<-make.W(dat,s.b)
#1st manipulation to get d(0)=1
n.c<-dim(H)[2]
n.r<-dim(H)[1]
P.c<-(P-H[,n.c])
H.c<-matrix(0,n.r,n.c-1)
for(i in 1:(n.c-1)){
H.c[,i]<-(H[,i]-H[,n.c])}
#2nd manipulation to get d(1/365)=rf
n.c2<-dim(H.c)[2]
n.r2<-dim(H.c)[1]
d.rf<-p(1,2,rf,0,0,0)/100
t<-(1/365)
k<-d.i
a<-a.i
P.c2<-(P.c-((d.rf-f(t,k,a))/(f(t,k-1,a)-f(t,k,a)))*H.c[,n.c2])
H.c2<-matrix(0,n.r2,n.c2-1)
for(i in 1:(n.c2-1)){
H.c2[,i]<-(H.c[,i]-((f(t,i,a)-f(t,k,a))/(f(t,k-1,a)-f(t,k,a)))*H.c[,n.c2])}
#get the coefficients back, price
Z.c2<-solve(t(H.c2)%*%W%*%H.c2)%*%t(H.c2)%*%W%*%P.c2
lkmo<-(d.rf-f(t,k,a))/(f(t,k-1,a)-f(t,k,a)) #The lambda_{k-1}
for(i in 1:length(Z.c2)){
lkmo<-(lkmo-Z.c2[i]*(f(t,i,a)-f(t,k,a))/(f(t,k-1,a)-f(t,k,a)))}
Z.c2<-c(Z.c2,lkmo)
Z.c2<-c(Z.c2,1-sum(Z.c2))
Phat<-H%*%Z.c2
Price.E<-round(P-Phat,2)
#yield
Yhat<-make.y(dates,Coupon,Phat,s.b,s.o,rf)
Yield.O<-Yield[2:length(Yield)]
Yield.E<-round(Yield.O-Yhat,3)
#Errors
mat<-matrix(0,2,2,dimnames = list(c("Cents", "Basis Points"), c("MAD", "RMSE")))
mat[1,1]<-sum(abs(100*Price.E))/length(Price.E)
mat[1,2]<-sqrt(sum((100*Price.E)^2)/length(Price.E))
mat[2,1]<-sum(abs(100*Yield.E))/length(Yield.E)
mat[2,2]<-sqrt(sum((100*Yield.E)^2)/length(Yield.E))
#plots
d.c<-d(Z.c2,l=l.i,a=a.i,j=j.i)
y.c<-yc(d.c,rf)
d2.c<-d2(Z.c2,l=l.i,a=a.i,j=j.i)
y2.c<-yc(d2.c,rf)
f.c<-fc(y.c,y2.c,rf)
#ouptup
out<-list(
coef=as.matrix(Z.c2),
price=data.frame(Phat,P,Price.E),
yield=data.frame(Yield.O,Yhat,Yield.E),
error=mat,
H=H,
W=W,
P=P,
H2=H.c2,
P2=P.c2,
alpha=a.i,
r=rf,
data=dat,
d=d.c,
y=y.c,
f=f.c
)
#output
out}

######################################################################
##########################run<-function###############################
######################################################################

run<-function(data,j.i=(1/12),l.i=10,OIS=TRUE, period=period){
#j.i=0.25;l.i=20;OIS=TRUE

RootDir <- getwd()
options(warn = -1)

#(a,b) the dimentions of the data matrix
a<-dim(data)[1]-1
b<-dim(data)[2]-1
data<-data[1:a,1:b]
nmat<-(l.i/j.i + 2)

#Matrix to store results.
BetaHat<-matrix(0,a-2,12,dimnames=list(as.character(data[3:a,1]),c("Beta1","Beta2","Beta3","Beta4","Beta5","Beta6","Beta7","Beta8","PriceErrorMAD", "YieldErrorMAD", "PriceErrorRMSE", "YieldErrorRMSE")))
DF<-matrix(0,a-2,nmat,dimnames=list(as.character(data[3:a,1]),c("Discount",round(seq(0,l.i,j.i),5))))
ZC<-matrix(0,a-2,nmat,dimnames=list(as.character(data[3:a,1]),c("Zero",round(seq(0,l.i,j.i),5))))
FR<-matrix(0,a-2,nmat,dimnames=list(as.character(data[3:a,1]),c("Forward",round(seq(0,l.i,j.i),5))))

#Construct the data matrix
for(q in 3:a){

#x is a vector of TRUE/FALSE saying if a bond is outstanding on that day
x<-!is.na(as.numeric(as.matrix(data[q,][2:b])))
#Number of outstanding bonds
n<-sum(x>0,na.rm=T)

Issue.ID<-rep(0,n)

#Coupon payments - only takes those applicable for the date
Coupon<-NULL
coup<-as.numeric(as.matrix(data[1,][2:b]))
j<-0
for(i in 1:b){
j<-j+1
ifelse(x[i],Coupon[j]<-coup[i],j<-j-1)}

#Maturity
ifelse(OIS==TRUE, {n1<-12; n2<-11; n3<-10}, {n1<-6; n2<-5; n3<-4})
mat<-as.Date(as.matrix(data[2,n1:b]),"%d/%m/%Y")#Maturities of bonds
matOIS<-as.numeric(as.matrix(data[2,2:n2]))	#Tenor in days of OIS
Maturity<-as.Date(data[q,1],"%d/%m/%Y")
#Sets the maturity of OIS
j<-0
for(i in 1:n3){
	j<-(j+1)
	ifelse(x[i],Maturity[j]<-Maturity[1]+matOIS[i],j<-(j-1))
}
#Sets maturity of bonds
j<-length(Maturity)
for(i in 1:b){
j<-(j+1)
ifelse(x[n3+i],Maturity[j]<-as.Date(mat[i]),j<-(j-1))
}

#Yields for relevant bonds
Yield<-NULL
yie<-as.numeric(as.matrix(data[q,][2:b]))
j<-0
for(i in 1:b){
j<-j+1
ifelse(x[i],Yield[j]<-yie[i],j<-j-1)}

#Make the data matrix
dat<-data.frame(Issue.ID,Coupon,Maturity,Yield)
rm(Issue.ID,Coupon,Maturity,Yield)

output<-go(dat,j.i=j.i,l.i=l.i)
BetaHat[q-2,1:8]=output$coef
BetaHat[q-2,9]<-output$error[1,1]
BetaHat[q-2,10]<-output$error[2,1]
BetaHat[q-2,11]<-output$error[1,2]
BetaHat[q-2,12]<-output$error[2,2]

DF[q-2,2:nmat]<-output$d[,2]
ZC[q-2,2:nmat]<-output$y[,2]
FR[q-2,2:nmat]<-output$f[,2]

print(q)}

setwd(paste(data_dir,"output",sep = "/"))
write.csv(BetaHat,paste(as.character(as.Date(data[3,1],"%d/%m/%Y")),"To",as.character(as.Date(data[a,1],"%d/%m/%Y")),"Beta",period,"csv",sep="."))
write.csv(round(DF,4),paste(as.character(as.Date(data[3,1],"%d/%m/%Y")),"To",as.character(as.Date(data[a,1],"%d/%m/%Y")),"Discount",period,"csv",sep="."))
write.csv(round(ZC,2),paste(as.character(as.Date(data[3,1],"%d/%m/%Y")),"To",as.character(as.Date(data[a,1],"%d/%m/%Y")),"Yield",period,"csv",sep="."))
write.csv(round(FR,2),paste(as.character(as.Date(data[3,1],"%d/%m/%Y")),"To",as.character(as.Date(data[a,1],"%d/%m/%Y")),"Forward",period,"csv",sep="."))

setwd(RootDir)
options(warn = 0)
}
