Dear all, I am using the following code for simulating Heston model. I am trying to compare skewness with respect to different Rhos, but it doesn't seem to work. Any body can tell me what is happening? I am using Euler-Maruyama Monte Carlo. Here is my sample code. Any hints would be really appreciated. theta=0.04 alpha=1.5 #delta=0.5 #rho=-0.9 n=100 m=1000 heston_sim<-function(rho,delta,m,mu){ s<-matrix(NA,n,m) v<-matrix(NA,n,m) for (j in 1:m){ for (i in 2:n) { s[1,j]=2 v[1,j]=0.4 dt=1/100 zv<-rnorm(1) zs<-rnorm(1) zx<-rho*zv+sqrt(1-rho^2)*zs w1<-zv*sqrt(dt) w2<-zx*sqrt(dt) s[i,j]=s[i-1,j]+(mu-0.5*max(v[i-1,j],0))*dt+sqrt(max(v[i-1,j],0))*w2 v[i,j]=v[i-1,j]+alpha*(theta-v[i-1,j])*dt+sqrt(max(v[i-1,j],0))*delta*w1 } } return(s) } s1<-heston_sim(-1,0.1,1000,0) s2<-heston_sim(0,0.1,1000,0) s3<-heston_sim(0.7,0.1,1000,0) a1<-c(diff(s1,lag=1)*100) a2<-c(diff(s2,lag=1)*100) a3<-c(diff(s3,lag=1)*100) plot(density(na.omit(a1),bw=0.5),main="",ylim=c(0,0.2)) lines(density(na.omit(a2),bw=0.5),col=2) lines(density(na.omit(a3),,bw=0.5),col=3) legend('topright',legend=c(expression(paste(rho, "=--0.3")),expression(paste(rho, "=0")),expression(paste(rho, "=0.5"))),col=1:3,lty=1) [[alternative HTML version deleted]]