n=1000;

vol = .25;
muoriginal = [-3.0642 -3.1160 -3.1886 -3.1009];
Foriginal = [0.047425 0.045743 0.043209 0.04791];
R = [0 0 0 0];
sigma = vol * (0.5)^(1/2);

R1 = zeros(1,n);
R2 = zeros(1,n);
R3 = zeros(1,n);
R4 = zeros(1,n);

for i=1:n
    mu = muoriginal; %reset mu and F
    F = Foriginal;
    
    for j=1:4
        R(j) = exp(randn*sigma + mu(j)); %generate R(j)
        e = log(R(j)) - log(F(j)); %calculate epsilon

        for k=j:4
            F(k) = exp(log(F(k)) + e);   %parallel shift
            mu(k) = log (F(k)) - sigma^2/2;   %recalc mu
        end

    end

    R1(i) = R(1); %store the R values
    R2(i) = R(2);
    R3(i) = R(3);
    R4(i) = R(4);

end

mean(R1) %display means and stds
mean(R2)
mean(R3)
mean(R4)
std(log(R1))
std(log(R2))
std(log(R3))
std(log(R4))