{

library(sm)
xdata=test[,6]
xdata=xdata*.0004690502

#Let us suppose that the data is in xdata
N=length(xdata)		#number of data points.

x1=min(xdata)
x2=max(xdata)
xeval=seq(x1,x2,length=100)
nevals=length(xeval)

#xeval is the array of 100 points at which you want to estimate
#the PDFs

zz=sm.density(xdata,model="Normal",eval.points=xeval)	#estimate the PDF
						#at the xeval points
xdensityorig=zz$estimate	#PDF of the original data

band=hnorm(xdata) / 2
band=hsj(xdata)

nsim=150

xsimpdf=matrix(0,nrow=nsim,ncol=nevals)

#generate bootstrap samples of the data and estimate the PDF
#of each bootstrap sample..

for(i in 1:nsim){

# you can replace this for monte carlo. E.g. x = rnorm(length(xdata), p1,p2)
# where p1 and p2 are the parameters of the normal distribution. Remember
# to comment out the boostrap command if using this..

#x=sample(xdata,replace=T)

#Smooth Bootstrap.

#x=sample(xdata,replace=T)+band*rnorm(N)

#varkern is the variance of the kernel you choose..
varkern=1		#variance of Normal kernel is 1

denom=sqrt(1 + (band*band*varkern/var(xdata)))

x=mean(xdata) + ((sample(xdata,replace=T) - mean(xdata) + band*rnorm(N))/denom)


# if you replace the above with the monte carlo then replace the
# following two command as well with 
# zz=dnorm(x,mean=mean(x),sd=sqrt(var(x)))
# xsimpdf[i,]=zz

zz=sm.density(x,eval.points=xeval,display="none")
xsimpdf[i,]=zz$estimate
}

# the following will producethe boxplot of the PDFs along
# with the PDF of the original data..

xs=1:nevals
xsimpdf1=as.data.frame(t(xsimpdf))
zz=boxplot(split(xsimpdf1,xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(xsimpdf,xdensityorig),xlab="",ylab="",cex=1.25)

z2=1:6
n1=1:6
z2[1]=z1[1]
z2[2]=z1[20]
z2[3]=z1[40]
z2[4]=z1[60]
z2[5]=z1[80]
z2[6]=z1[100]

n1[1]=xeval[1]
n1[2]=xeval[20]
n1[3]=xeval[40]
n1[4]=xeval[60]
n1[5]=xeval[80]
n1[6]=xeval[100]
n1=round(n1,dig=0)
n1=as.character(n1)

axis(1,at=z2,labels=n1,cex=1.00)
lines(z1,xdensityorig,lty=2,lwd=2)

title(main="PDFs from the simulations and the historical data")
}