Monte Carlo simulation with "double" indicization

Good evening everybody,

I would like run a Monte Carlo simulation with the following code. My problem is that I have to define kind of "double" loop, one for the time "t" and the other one for the "n" realization of the random variable rn. In the following my code WITHOUT Monte Carlo sampling (it works fine).

T=10000;
sigma=0.157;
k=0.3;
beta=0.99;
po=0.35;
phipi=1.5;
phiz=0;   
B2=[sigma 1-beta*phipi; k*sigma k+beta*(sigma+phiz)]  
B1=1/(sigma+phiz+k*phipi)                           
B=B1*B2

y=NaN(2,T)     

y(1,1)=0.1    
y(2,1)=0.2    

kappa=[1/(sigma+phiz+phipi*k); k/(sigma+phiz+phipi*k)]    
rn=randn(1,T)

a=NaN(2,T);
a(1,1)=1
a(2,1)=2


for t=2:T
y(:,t)=B*a(:,t-1)+kappa*rn(:,t);                          
a(:,t)=a(:,t-1)+t^(-1)*(y(:,t)-a(:,t-1));             
end

stz=a(1,:)'             
stpi=a(2,:)'  

In the following instead "what I am trying to do". Bear in mind that I would like to obtain a multidimensional vector defined as y(2,T,N) in order to plot a "density function" of the realizations of the components of vector y (that is, the two variables z e pi).

T=100;
sigma=0.157;
k=0.3;
beta=0.99;
po=0.35;
phipi=1.5;
phiz=0;   
B2=[sigma 1-beta*phipi; k*sigma k+beta*(sigma+phiz)]  
B1=1/(sigma+phiz+k*phipi)                           
B=B1*B2

y=NaN(2,T)     

y(1,1)=0.1    
y(2,1)=0.2    

kappa=[1/(sigma+phiz+phipi*k); k/(sigma+phiz+phipi*k)]    
rn=randn(1,T)

a=NaN(2,T);
a(1,1)=1
a(2,1)=2

N=100
for n=1:N 
    rn=randn(1,N)
    for t=2:T
        y(:,t,n)=B*a(:,t-1)+kappa*rn(:,t);                           
        a(:,t)=a(:,t-1)+t^(-1)*(y(:,t)-a(:,t-1));                
    end
end

Am I doing well? Another problem is that I cannot push the simulation above N=100 or T=100 (it is seems to me that N e T must agree, am I right?) otherwise it take too much time to compute.

What command I am supposed to use in Matlab to plot a density function of the MC simulation?

Thank you everybody

[Image Analyst]

Not sure I understand the question, and I don't see any variable "m". Are you trying to run an experiment some number of times (how many?), where t takes on a random value between 0 and T, and n takes on a random value between 1 and N, and your two equations for y and a depend on those random values?

Why do you use y as a 3-D variable when you assign y, but as a 2-D variable when you use y to compute "a"?

[Uccio]

Thank you for your response,

rn is just defined as "rn=randn(1,N)". You can find it written inside the only loop.

I am trying to run the experiment N=100 times. t does not take random value, it is simple the time, that is, i am saying to matlap to repeat the loop for t=2, t=3, t=4 and so on. The same holds for n: I am saying matlap to repeat the loop for the first realziation of rn, for the second one etc..

I quote you "Why do you use y as a 3-D variable when you assign y, but as a 2-D variable when you use y to compute "a"?". THAT is exactly the object of my question: I am having hard time in writing a montecarlo simulation and I would like to know if my point is correct or not.

Thank you again for help

[John D'Errico]

A bit confusing. You say this:

"In the following my code WITHOUT Monte Carlo sampling (it works fine)."

Yet then you have the line:

rn=randn(1,T)

Which makes this indeed a Monte Carlo sampling method.

By the way, learn to use semi-colons! If not, then you will dump out 10000 numbers to the screen every time you run this code.

[Uccio]

I get your point. Actually rn(t) cannot be a scalar (I cannot write rn=randn(1)) since rn(t) is a variable in a linear dynamic sthocastic model. At each point in time, I allow it to have a different realization.

Now my question is: cannot I run my model N times and then aggregate the results in order to obtain a density function for z and pi?