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: [email protected] > Sent: Thu, 4 Apr 2013 15:14:05 +0800 > To: [email protected] > 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]] > > ______________________________________________ > [email protected] 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! ______________________________________________ [email protected] 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.
