require(fOptions)
require(MASS)
multiassetpath<-function(S0,drift,TVal,covmat,NSteps,NRepl=1,NAssets)
{
#S0 is spot array of length NAssets
# drift is drift array of length NAssets
# covmat is the variance covariance matrix (NAssets x NAssets)
# NSteps is number of time steps (needed if your payout is path dependent)
# NRepl is number of simulations 
# Nassets is the number of assets.
# returns an array of size NSteps x NPaths x NAsset
Path<-array(0,c(NSteps,NRepl,NAssets))
DeltaT<-TVal/NSteps

# some alternative for generating random numbers
# randmatrix with pseudo random numbers 
#randmat <- matrix(rnorm(NAssets*NRepl*NSteps),ncol=NAssets) %*% chol(covmat)
#randmat<-mvrnorm(NRepl*NSteps,c(0,0),Sigma=covmat)
# using quasi-random numbers sobol and halton
#randmat <- matrix(rnorm.halton(NAssets*NRepl*NSteps,NAssets),ncol=NAssets) %*% chol(covmat)
#randmat <-matrix(rnorm.sobol(NRepl,NAssets*NSteps),ncol=NAssets*NSteps) 
Shocks<-array(0,c(NSteps+1,NRepl,NAssets))
for(ii in 1:NAssets)
{
Shocks[1,,ii]<-S0[ii]
}
mudt<-drift*DeltaT
randmat <-matrix(rnorm.sobol(NRepl,NAssets*NSteps),ncol=NAssets*NSteps) 
cholmat<-chol(covmat)
for(kk in 1:NRepl)
{
#sobol
myrandmat<-matrix(randmat[kk,],nrow=NSteps,ncol=NAssets) %*% cholmat

for(zz in 1:NSteps)
{
for(jj in 1:NAssets) {
Shocks[zz+1,kk,jj]<-Shocks[zz,kk,jj]*(exp(mudt[jj]+myrandmat[zz,jj]))
Path[zz,kk,jj]<-Shocks[zz+1,kk,jj]
}
}
}

return(Path)
}

