此处 t 代表行权时间
r代表连续复利收益率    sig 波动率      m代表模拟m次      n代表每次模拟n步
div代表连续股利收益率
我下面提供的是部分欧式看涨期权的蒙特卡洛模拟的核心代码  如有问题请纠正
但是我本人sas初学 但知道其有强大的分析和作图功能
请高手将我的代码丰富下 能够实现图像形式的股价模拟路径 非常感谢



dt = t/n;                                                                                                                           
nudt = (r-div-0.5*sig^2)*dt;                                                                                                         
sigsdt = sig*sqrt(dt);                                                                                                               
lns = ln(s);                                                                                                                                                                                                                                      
sum_CT = 0;                                                                                                                          
sum_CT2 = 0;                                                                                                                           
do  j = 1 to M;     *for each simulation;                                                                                            
lnst = lns;                                                                                                                           
  do i = 1 to  N ;     *for each time step;                                                                                       
  X=rannor(0);                                                                                                                     
  lnst=lnst+nudt+sigsdt*x;   *evolve the stock price;                                                                              
  end;                                                                                                         

ST = exp(lnst);                                                                                                                     
CT = max{0,ST-K};                                                                                                                    
sum_CT = sum_CT + CT;                                                                                                                 
sum_CT2 = sum_CT2 + CT*CT;                                                                                                            
end;     

call_value = sum_CT/M*exp(-r*t);                                                                                                  
SD = sqrt((sum_CT2 - sum_CT*sum_CT/M)*exp(-2*r*T)/(M-1));                                                                             
SE = sd/sqrt(M);