Thank you, Sander. Exp outside and Random Number generation with specified mean and standard deviation is a nice tip.

I don't know very well how to write Mathematica code to be compiled in C. Nevertheless, I tried that below, but it does not improve much. What am I doing wrong here?

tryC = Compile[{{n, _Integer}},
   Block[{spot = ConstantArray[95.0, n], K = 100, \[Sigma] = 0.15, 
     T = 1, r = 0.01, H = 85, m = 50, h, \[Alpha], RND, index, cc},
    h = T/m;
    \[Alpha] = r - \[Sigma]^2/2;
    index = ConstantArray[0.0, n];
    RND = 
     RandomVariate[
      NormalDistribution[\[Alpha] h, \[Sigma] Sqrt[h]], {m, n}];
    Table[spot = spot RND[[i]];
     cc = Clip[spot, {H + 33 $MachineEpsilon, H}, {i, 0}];
     index += Clip[index + cc, {i, i}, {0, 0}];
     , {i, m}];
    index],
   CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
   "CompilationTarget" -> "C", "RuntimeOptions" -> "Speed"];

index = IntegerPart@
    Clip[tryC[1000000], {1, 50}, {50, 50}]; // AbsoluteTiming

Another question is why does it take so long to generate the random numbers? Is is because these numbers are "more random"? or because they are more precise (more decimal digits)?

In your updated code, if you do the scaling and Exp outside, it speeds up like 40% or so on my laptop:

spot=S;index=ConstantArray[0,n]; (* vector specifying the hitting times for each path *)
RND=Exp[RandomVariate[NormalDistribution[\[Alpha] h,\[Sigma] Sqrt[h]],{m,n}]];
Table[spot*= RND[[i]]; (* simulation of prices in period i *)
    cc=Clip[spot,{H+33$MachineEpsilon,H},{i,0}];
    index+=Clip[index+cc,{i,i},{0,0}]; (* record the hitting times to index *)
,{i,m}];//AbsoluteTiming
index=IntegerPart@Clip[index,{1,50},{50,50}];//AbsoluteTiming

In my matlab 'converted' code, everything is a PackedArray (Developer`PackedArrayQ), so also there, there is not much room for speedups...

I'm not sure if you can really make it a lot faster in Mathematica to be honest without resorting to C and/or parallelisation.

You can translate your Matlab code to:

   S0=95.0;
   K=100;
   sigma=0.15;
   T=1;
   r=0.01
   alpha=(r-1/2*sigma^2);
   H=85;
   b=-Log[H/S0];
   c=Log[K/S0];
   m=50;
   h=T/m;

   n=1000000;
   AbsoluteTiming[mZ =Exp[RandomVariate[NormalDistribution[alpha h,sigma Sqrt[h]],{m,n}]];]

   vS=ConstantArray[S0,n];
   vtau=ConstantArray[m,n];

   Do[
       vS *= mZ[[i]];
       tmp=UnitStep[H-vS](1-Unitize[vtau-m]);
       vtau*= tmp;
       vtau+= tmp i;
   ,
       {i,m}
   ]//AbsoluteTiming

Most time is in the creation of the random numbers, which you probably can't really speed up...

Not always, Sander. If you have a list with some of the elements equal to i, then it will give zero everywhere except those elements which equal i.

For instance, let i = 3, then

Clip[{0, 1, 3, 6, 3, 9, 3}, {3, 3}, {0, 0}]

gives:

{0, 0, 3, 0, 3, 0, 3}

I figured this is a faster way to do a Replace operation.

So, in the code, this line:

cc = Clip[spot, {H + 33 $MachineEpsilon, H}, {i, 0}];

identifies which which paths go below H and gives a vector cc with time i in positions corresponding to these paths and zero for others.

Then, this line of code:

index += Clip[index + cc, {i, i}, {0, 0}];

adds to index only the paths which hit the barrier in time i, eliminating those which did that previously ((a previous hitting time + i) > i replaced by zero in the resulting Clip vector).

I see now that the price process simulation alone takes more than the entire MATLAB code. Why are you saying that it is not vectorized?