library("quadprog")

##############################Preparing for data
rawdata <- read.csv ("C:/Adjusted US returns and Adjusted t-bills.csv",header=T)
stocks = rawdata[,2]
bonds = rawdata[,1]
data = data.frame(bonds,stocks)
OriData = as.matrix(data)

##############################the GetBSData function
GetBSData<-function(data){
x = 1:396
s = sample(x,6,replace=T)
bsdata = data[(s[1]):(s[1]+59),]
	for (j in 2:6) {
  		a = data[(s[j]):(s[j]+59),]
  		bsdata = rbind(bsdata,a)
	}
return(bsdata)
}

#set.seed(1234)
#trial<-GetBSData(OriData)

##############################the Minimisation function
Opt<-function(data, horizon, lamda){
StockReturn<-numeric(30/horizon)
BondReturn<-numeric(30/horizon)
for (x in 1: (30/horizon)){
	StockReturn[x]<-prod(data[(12*horizon*(x-1)+1):(12*horizon*(x-1)+12*horizon),2])-1
	BondReturn[x]<-prod(data[(12*horizon*(x-1)+1):(12*horizon*(x-1)+12*horizon),1])-1
}
Return<-cbind(BondReturn,StockReturn)
MeanVec<-c(mean(BondReturn),mean(StockReturn))
VCovMat<-cov(Return)
#return(MeanVec, VCovMat)

a<-c(1,1)
a<-cbind(a, diag(1,2))

WtVec<-solve.QP(Dmat=VCovMat*2, dvec= MeanVec*lamda,Amat=a,bvec=c(1,0,0),meq=1)

#return(MeanVec, VCovMat, WtVec$solution)
return(WtVec$solution)
}

#Opt(OriData+1, 5, 0)
Opt(OriData+1, 5, 0)

##############################
set.seed(4114)
bs=1000						###number of bootstrap samples
lamdaseq<-seq(0,1,.05)				###the lamda sequence. currently from 0 to 1 by .05.

x<-numeric(bs*length(lamdaseq))		###
w1<-matrix(x, bs, length(lamdaseq))		###To initialise the matrices.
w5<-matrix(x, bs, length(lamdaseq))		###1, 5, 10 denote the horizon.
w10<-matrix(x, bs, length(lamdaseq))	###

for (i in 1: bs){
BSData<-GetBSData(OriData)+1
j=1
	for (lamda in lamdaseq){
		w1[i,j]<-Opt(BSData, 1, lamda)[1]
		w5[i,j]<-Opt(BSData, 5, lamda)[1]
		w10[i,j]<-Opt(BSData, 10, lamda)[1]
		j=j+1
	}

x<-numeric(length(lamdaseq)*9)		###To initialise the table
table<-matrix(x, length(lamdaseq), 9)	###

for (k in 1:length(lamdaseq)){		#k:index for lamda

table[k,1]<-sort(w1[,k])[.05*bs]		###The first 3 cols are for 1-yr horizon.
table[k,2]<-mean(w1[,k])			###From left to right: 5 percentile,
table[k,3]<-sort(w1[,k])[.95*bs]		###mean, and 95 percentile.

table[k,4]<-sort(w5[,k])[.05*bs]		###
table[k,5]<-mean(w5[,k])			###Col 4-6 are for 5-yr horizon.
table[k,6]<-sort(w5[,k])[.95*bs]		###

table[k,7]<-sort(w10[,k])[.05*bs]		###
table[k,8]<-mean(w10[,k])			###Col 7-9 are for 5-yr horizon.
table[k,9]<-sort(w10[,k])[.95*bs]		###
}
}

table

TenMinusOne<-numeric(length(lamdaseq))
FiveMinusOne<-numeric(length(lamdaseq))
TenMinusFive<-numeric(length(lamdaseq))

for (p in 1:length(lamdaseq)){
DiffVec<-w10[,p]-w1[,p]
TenMinusOne[p]<-length(DiffVec[DiffVec>0])

DiffVec<-w5[,p]-w1[,p]
FiveMinusOne[p]<-length(DiffVec[DiffVec>0])

DiffVec<-w10[,p]-w5[,p]
TenMinusFive[p]<-length(DiffVec[DiffVec>0])
}

diff<-cbind(FiveMinusOne,TenMinusOne)
diff<-cbind(diff, TenMinusFive)
sn<-seq(1, length(lamdaseq))
f2<-cbind(sn, diff)
f2

##############################################END