As attached is the txt file for the code.
Regards,BL
> Date: Thu, 4 Apr 2013 03:19:47 -0800
> From: jrkrid...@inbox.com
> Subject: RE: [R] Help for bootstrapping
> To: boon_lo...@hotmail.com; r-help@r-project.org
>
> It looks like you formatted the code in html and it is essentially impossible
> to read. Can you resend in plain text?
>
> Thanks
>
> John Kane
> Kingston ON Canada
>
>
> > -----Original Message-----
> > From: boon_lo...@hotmail.com
> > Sent: Thu, 4 Apr 2013 15:14:05 +0800
> > To: r-help@r-project.org
> > Subject: [R] Help for bootstrapping
> >
> > I have a set of data for US t-bill returns and US stock returns frm
> > 1980-2012. I am trying to bootstrap the data and obtain the minimum
> > variance portfolio and repeat this portfolio 1000 times. However I am
> > unable to get the correct code function for the minimum variance
> > portfolio. When I tried to enter Opt(OriData+1, 1, 5, 0), I get
> > "error:subscript out of bounds" Please help!
> > library("quadprog")
> > ##############################Preparing for datarawdata =
> > read.table("C:/Desktop/data.txt", header=T)Rf = rawdata[,1]US =
> > rawdata[,2]data = data.frame(Rf,US)OriData = as.matrix(data)
> > ##############################the GetBSData
> > functionGetBSData<-function(data){x = 1:396s =
> > 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
> > functionOpt<-function(data, horizon, col,
> > lamda){TbillReturn<-numeric(30/horizon)USReturn<-numeric(30/horizon)for
> > (x in 1: (30/horizon)){
> TbillReturn[x]<-prod(data[(12*horizon*(x-1)+1):(12*horizon*(x-1)+12*horizon),col])-1
> USReturn[x]<-prod(data[(12*horizon*(x-1)+1):(12*horizon*(x-1)+12*horizon),2])-1}Return<-cbind(TbillReturn,USReturn)MeanVec<-c(mean(TbillReturn),mean(USReturn))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, 1, 5, 0)
> > ##############################set.seed(4114)bs=1000
> > ###number of
> > bootstrap samplesRegion<-5
> > ###Region indecies, check
> > above.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)+1j=1 for (lamda in lamdaseq){
> > w1[i,j]<-Opt(BSData, 1, Region, lamda)[1] w5[i,j]<-Opt(BSData, 5,
> > Region, lamda)[1] w10[i,j]<-Opt(BSData, 10, Region, lamda)[1]
> > j=j+1 }
> > x<-numeric(length(lamdaseq)*9) ###To initialise the
> > tabletable<-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
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help@r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
> ____________________________________________________________
> FREE 3D EARTH SCREENSAVER - Watch the Earth right on your desktop!
> Check it out at http://www.inbox.com/earth
>
>
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______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
