{

#simulate from Empirical CDF - i.e. BOOTSTRAP 
# or from a Kernel density PDF - i.e., SMOOTHED BOOTSTRAP

data=matrix(scan("aismr.txt"),ncol=13,byrow=T)
xdata=data[,11]

#suppose the historical data is in the vector 'xdata'
N=length(xdata)		#number of data points.

#Monte Carlo...
nsim=150	#number of simulations, each of length N

#points at which to estimate the PDF from the simulations and the observed

xeval=seq(min(xdata)-sd(xdata),max(xdata)+sd(xdata),length=100)
neval=length(xeval)
xsimpdf=matrix(0,nrow=nsim,ncol=neval)

#initialize vectors for the statistics
meansim=1:nsim
mediansim=1:nsim
sdsim=1:nsim
skewsim=1:nsim
iqrsim=1:nsim

# Evaluate the Kernel density PDF at the evaluation points based on the data
library(sm)
xdensityorig=sm.density(xdata, eval.points=xeval, display="none")$estimate

#compute bandwidth of the data
band=hnorm(xdata)
######## Simulation ##########


for(i in 1:nsim){

# Simulate from the desired model (i.e., PDF) 
# each simulation should of the same length as the original data
# used to fit the model, N = length(xdata)

#simulate from empirical CDF - BOOTSTRAP
N = length(xdata)
xsim=sample(xdata,N, replace=TRUE)

#From Kernel PDF  - smooth bootstrap..
#xsim=sample(xdata,N, replace=TRUE) + band*rnorm(N)

#From Kernel PDF - Smooth bootstrap with correction to reproduce the variance

#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)))
#xsim=mean(xdata) + ((sample(xdata,N,replace=TRUE) - mean(xdata) + band*rnorm(N))/denom)

# compute the statistics from the simulation
meansim[i]=mean(xsim)
sdsim[i]=sd(xsim)
#sdsim[i]=sd(log(xsim))
skewsim[i]=skew(xsim)
iqrsim[i]=diff(quantile(xsim,c(0.25,0.75)))
mediansim[i]=quantile(xsim,c(0.5))

#estimate the Kernel PDF at the evaluation points based on the simulated data

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

}


#boxplot the stats from the simulations and overlay the corresponding value
# from the original data as a point.

par(mfrow=c(1,5))

zz=boxplot(meansim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,mean(xdata),col="red",cex=1.25,pch=16)
title(main="Mean")

zz=boxplot(mediansim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,quantile(xdata,c(0.5)),col="red",cex=1.25, pch=16)
title(main="Median")

zz=boxplot(sdsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
#points(z1,sd(log(xdata)),col="red",cex=1.25)
points(z1,sd(xdata),col="red",cex=1.25,pch=16)
title(main="SD")

zz=boxplot(skewsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,skew(xdata),col="red",cex=1.25,pch=16)
title(main="Skew")

zz=boxplot(iqrsim,xlab="",ylab="Mean",plot=F)
zz$names=rep("",length(zz$names))
z1=bxp(zz,xlab="",ylab="Mean",cex=1.25)
points(z1,diff(quantile(xdata,c(0.25,0.75))),col="red",cex=1.25,pch=16)
title(main="IQR")

par(ask=TRUE)

#boxplot the PDFs from the simulations along with that of the original data

par(mfrow=c(1,1))
xs=1:neval

# For R version  2.4.1 
zz=boxplot(split(t(xsimpdf),xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(xsimpdf,xdensityorig),xlab="",ylab="",cex=1.25)

npts=10         #number of points to plot on the x-axis..
n2=round(neval*(1:npts)/npts)

z2=z1[n2]
n1=xeval[n2]


n1=round(n1,dig=2)
n1=as.character(n1)

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

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


}