Multi-Thread calling a function

The reason for the crash is that the mat_vec procedure which is called at the end of the threadfor loop is assigning to many different global variables.

We can help you to make this work, and most likely simplify this code as well. However, to do that I have some questions about the code.

I don't see any use of the variables which are assigned to in the mat_vec procedure. For example, I see these assignments:

hp2_1pt_2pt= arraytomat(hp2_1pt_2pt);
hp2_1ft_2pt= arraytomat(hp2_1ft_2pt);
hp2_1np_2ft= arraytomat(hp2_1np_2ft);
hp2_1ue_2ft= arraytomat(hp2_1ue_2ft);
hp2_1pt_2ft= arraytomat(hp2_1pt_2ft);
hp2_1ft_2ft= arraytomat(hp2_1ft_2ft);

However, I don't see them referenced in or after the threadfor loop. Do you need them during the loop, or do you just need a print out of the means, etc AFTER the threadfor loop? ...or maybe something else I did not think of. 🙂

Well actually I need them during the procedure. because I will use them to compare decisions. I will show you.

All of these is inside a procedure:

// Hosehold's types (observed characteristics)
//Other member & Kids & Non Income
//1 NO NO NO
//2 NO NO YES
//3 NO YES NO
//4 NO YES YES
//5 YES NO NO
//6 YES NO YES
//7 YES YES NO
//8 YES YES YES
 
// 1=Have children 2=Not have children
kids_matrix1=1*((HHtype_vect.==3)+(HHtype_vect.==4)+(HHtype_vect.==7)+(HHtype_vect.==8))+ 2*((HHtype_vect.==1)+(HHtype_vect.==2)+(HHtype_vect.==5)+(HHtype_vect.==6));
// 1=Other HH members>12 2=No-other HH members>12
d13_matrix= 1*((HHtype_vect.==5)+(HHtype_vect.==6)+(HHtype_vect.==7)+(HHtype_vect.==8))+ 2*((HHtype_vect.==1)+(HHtype_vect.==2)+(HHtype_vect.==3)+(HHtype_vect.==4));
if CC_time==1;
    // 1=Other HH members>12 and No-other HH members+Children 2=No-other HH members>12
    d13_matrix= 1*((HHtype_vect.==5)+(HHtype_vect.==6)+(HHtype_vect.==7)+(HHtype_vect.==8)+(HHtype_vect.==3)+(HHtype_vect.==4))+ 2*((HHtype_vect.==1)+(HHtype_vect.==2));
endif;
// Discrete binary distribution
s1 = sea_cost ;
s2 = sea_cost ;
simul_s1 = 1*(zeta_s_m.<=gam_s1) + 2*(zeta_s_m.>gam_s1) ; // gam_s1 is the prob of low searching cost
simul_s2 = 1*(zeta_s_f.<=gam_s2) + 2*(zeta_s_f.>gam_s2) ; // gam_s1 is the prob of low searching cost
meanc(simul_s1)~meanc(simul_s2);
//productivity vector
prodM_vector = 1 | prod_h ;
prodW_vector = 1 | prod_h ;
// Wage Distribution
PDF_ln = Wage_PDF(1);
// CDF of wage offer distribution and SE distribution
P_M_pt_cum = cumsumc(PDF_ln.w1_pt*step_length_w1);
P_M_pt_cum = P_M_pt_cum./P_M_pt_cum[N_w];//=1; // CDF of men's wage grid
P_M_ft_cum = cumsumc(PDF_ln.w1_ft*step_length_w1);
P_M_ft_cum = P_M_ft_cum./ P_M_ft_cum[N_w];//=1; // CDF of men's wage grid
P_W_pt_cum = cumsumc(PDF_ln.w2_pt*step_length_w1);
P_W_pt_cum = P_W_pt_cum./ P_W_pt_cum[N_w];//=1; //CDF of women's wage grid
P_W_ft_cum = cumsumc(PDF_ln.w2_ft*step_length_w1);
P_W_ft_cum = P_W_ft_cum./P_W_ft_cum[N_w];//=1; //CDF of women's wage grid
//P_M_pt_cum~P_M_ft_cum~P_W_pt_cum~P_W_ft_cum;
 
// Value Function Fixed Point for each Non-Labor Income hasta aqui vamos 11/26
for dHH (1,2,1);
    for k (1,2,1); //males
        for l (1,2,1); //females
            time_in_s=hsec;
            //s1 = s1_vector[k] ;
            //s2 = s2_vector[l] ;
            prod1=prodM_vector[k];
            prod2=prodW_vector[l];
            //print "Search Cost: s1" k "s2" l "prod1" prod1 "prod2" prod2 "d13" dHH ;
            // Non-Labor Income 1
            time_in=hsec;
            V.V_1np_2np = arraytomat(V_nl.V_1np_2np[1,dHH,k,l,.,.]) ;
            V.V_1ue_2np = arraytomat(V_nl.V_1ue_2np[1,dHH,k,l,.,.]) ;
            V.V_1np_2ue = arraytomat(V_nl.V_1np_2ue[1,dHH,k,l,.,.]) ;
            V.V_1ue_2ue = arraytomat(V_nl.V_1ue_2ue[1,dHH,k,l,.,.]) ;
            V.V_1np_2pt = arraytomat(V_nl.V_1np_2pt[1,dHH,k,l,.,.]) ;
            V.V_1ue_2pt = arraytomat(V_nl.V_1ue_2pt[1,dHH,k,l,.,.]) ;
            V.V_1np_2ft = arraytomat(V_nl.V_1np_2ft[1,dHH,k,l,.,.]) ;
            V.V_1ue_2ft = arraytomat(V_nl.V_1ue_2ft[1,dHH,k,l,.,.]) ;
 
            V.V_1pt_2np = arraytomat(V_nl.V_1pt_2np[1,dHH,k,l,.,.]) ;
            V.V_1ft_2np = arraytomat(V_nl.V_1ft_2np[1,dHH,k,l,.,.]) ;
            V.V_1pt_2ue = arraytomat(V_nl.V_1pt_2ue[1,dHH,k,l,.,.]) ;
            V.V_1ft_2ue = arraytomat(V_nl.V_1ft_2ue[1,dHH,k,l,.,.]) ;
            V.V_1pt_2pt = arraytomat(V_nl.V_1pt_2pt[1,dHH,k,l,.,.]) ;
            V.V_1ft_2pt = arraytomat(V_nl.V_1ft_2pt[1,dHH,k,l,.,.]) ;
            V.V_1pt_2ft = arraytomat(V_nl.V_1pt_2ft[1,dHH,k,l,.,.]) ;
            V.V_1ft_2ft = arraytomat(V_nl.V_1ft_2ft[1,dHH,k,l,.,.]) ;
 
            V.V_1np_2np_nk = arraytomat(V_nl.V_1np_2np_nk[1,dHH,k,l,.,.]);
            V.V_1ue_2np_nk = arraytomat(V_nl.V_1ue_2np_nk[1,dHH,k,l,.,.]);
            V.V_1np_2ue_nk = arraytomat(V_nl.V_1np_2ue_nk[1,dHH,k,l,.,.]);
            V.V_1ue_2ue_nk = arraytomat(V_nl.V_1ue_2ue_nk[1,dHH,k,l,.,.]);
            V.V_1np_2pt_nk = arraytomat(V_nl.V_1np_2pt_nk[1,dHH,k,l,.,.]);
            V.V_1ue_2pt_nk = arraytomat(V_nl.V_1ue_2pt_nk[1,dHH,k,l,.,.]);
            V.V_1np_2ft_nk = arraytomat(V_nl.V_1np_2ft_nk[1,dHH,k,l,.,.]);
            V.V_1ue_2ft_nk = arraytomat(V_nl.V_1ue_2ft_nk[1,dHH,k,l,.,.]);
 
            V.V_1pt_2np_nk = arraytomat(V_nl.V_1pt_2np_nk[1,dHH,k,l,.,.]);
            V.V_1ft_2np_nk = arraytomat(V_nl.V_1ft_2np_nk[1,dHH,k,l,.,.]);
            V.V_1pt_2ue_nk = arraytomat(V_nl.V_1pt_2ue_nk[1,dHH,k,l,.,.]);
            V.V_1ft_2ue_nk = arraytomat(V_nl.V_1ft_2ue_nk[1,dHH,k,l,.,.]);
            V.V_1pt_2pt_nk = arraytomat(V_nl.V_1pt_2pt_nk[1,dHH,k,l,.,.]);
            V.V_1ft_2pt_nk = arraytomat(V_nl.V_1ft_2pt_nk[1,dHH,k,l,.,.]);
            V.V_1pt_2ft_nk = arraytomat(V_nl.V_1pt_2ft_nk[1,dHH,k,l,.,.]);
            V.V_1ft_2ft_nk = arraytomat(V_nl.V_1ft_2ft_nk[1,dHH,k,l,.,.]);
 
            time_in=hsec;
            U = hh_utility(nonlabInc1,dHH);
 
            //HOUSEHOLD PRODUCTION MALE HOURS
            hp1.V_1np_2np[1,dHH,k,l,.,.] = U.hh1_1np_2np ;
            hp1.V_1np_2ue[1,dHH,k,l,.,.] = U.hh1_1np_2ue ;
            hp1.V_1np_2pt[1,dHH,k,l,.,.] = U.hh1_1np_2pt ;
            hp1.V_1np_2ft[1,dHH,k,l,.,.] = U.hh1_1np_2ft ;
            hp1.V_1ue_2np[1,dHH,k,l,.,.] = U.hh1_1ue_2np ;
            hp1.V_1ue_2ue[1,dHH,k,l,.,.] = U.hh1_1ue_2ue ;
            hp1.V_1ue_2pt[1,dHH,k,l,.,.] = U.hh1_1ue_2pt ;
            hp1.V_1ue_2ft[1,dHH,k,l,.,.] = U.hh1_1ue_2ft ;
            hp1.V_1pt_2np[1,dHH,k,l,.,.] = U.hh1_1pt_2np ;
            hp1.V_1pt_2ue[1,dHH,k,l,.,.] = U.hh1_1pt_2ue ;
            hp1.V_1pt_2pt[1,dHH,k,l,.,.] = U.hh1_1pt_2pt ;
            hp1.V_1pt_2ft[1,dHH,k,l,.,.] = U.hh1_1pt_2ft ;
            hp1.V_1ft_2np[1,dHH,k,l,.,.] = U.hh1_1ft_2np ;
            hp1.V_1ft_2ue[1,dHH,k,l,.,.] = U.hh1_1ft_2ue ;
            hp1.V_1ft_2pt[1,dHH,k,l,.,.] = U.hh1_1ft_2pt ;
            hp1.V_1ft_2ft[1,dHH,k,l,.,.] = U.hh1_1ft_2ft ;
 
            hp1.V_1np_2np_nk[1,dHH,k,l,.,.] = U.hh1_1np_2np_nk;
            hp1.V_1np_2ue_nk[1,dHH,k,l,.,.] = U.hh1_1np_2ue_nk;
            hp1.V_1np_2pt_nk[1,dHH,k,l,.,.] = U.hh1_1np_2pt_nk;
            hp1.V_1np_2ft_nk[1,dHH,k,l,.,.] = U.hh1_1np_2ft_nk;
            hp1.V_1ue_2np_nk[1,dHH,k,l,.,.] = U.hh1_1ue_2np_nk;
            hp1.V_1ue_2ue_nk[1,dHH,k,l,.,.] = U.hh1_1ue_2ue_nk;
            hp1.V_1ue_2pt_nk[1,dHH,k,l,.,.] = U.hh1_1ue_2pt_nk;
            hp1.V_1ue_2ft_nk[1,dHH,k,l,.,.] = U.hh1_1ue_2ft_nk;
            hp1.V_1pt_2np_nk[1,dHH,k,l,.,.] = U.hh1_1pt_2np_nk;
            hp1.V_1pt_2ue_nk[1,dHH,k,l,.,.] = U.hh1_1pt_2ue_nk;
            hp1.V_1pt_2pt_nk[1,dHH,k,l,.,.] = U.hh1_1pt_2pt_nk;
            hp1.V_1pt_2ft_nk[1,dHH,k,l,.,.] = U.hh1_1pt_2ft_nk;
            hp1.V_1ft_2np_nk[1,dHH,k,l,.,.] = U.hh1_1ft_2np_nk;
            hp1.V_1ft_2ue_nk[1,dHH,k,l,.,.] = U.hh1_1ft_2ue_nk;
            hp1.V_1ft_2pt_nk[1,dHH,k,l,.,.] = U.hh1_1ft_2pt_nk;
            hp1.V_1ft_2ft_nk[1,dHH,k,l,.,.] = U.hh1_1ft_2ft_nk;
 
            //HOUSEHOLD PRODUCTION FEMALE HOURS
            hp2.V_1np_2np[1,dHH,k,l,.,.] = U.hh2_1np_2np ;
            hp2.V_1np_2ue[1,dHH,k,l,.,.] = U.hh2_1np_2ue ;
            hp2.V_1np_2pt[1,dHH,k,l,.,.] = U.hh2_1np_2pt ;
            hp2.V_1np_2ft[1,dHH,k,l,.,.] = U.hh2_1np_2ft ;
            hp2.V_1ue_2np[1,dHH,k,l,.,.] = U.hh2_1ue_2np ;
            hp2.V_1ue_2ue[1,dHH,k,l,.,.] = U.hh2_1ue_2ue ;
            hp2.V_1ue_2pt[1,dHH,k,l,.,.] = U.hh2_1ue_2pt ;
            hp2.V_1ue_2ft[1,dHH,k,l,.,.] = U.hh2_1ue_2ft ;
            hp2.V_1pt_2np[1,dHH,k,l,.,.] = U.hh2_1pt_2np ;
            hp2.V_1pt_2ue[1,dHH,k,l,.,.] = U.hh2_1pt_2ue ;
            hp2.V_1pt_2pt[1,dHH,k,l,.,.] = U.hh2_1pt_2pt ;
            hp2.V_1pt_2ft[1,dHH,k,l,.,.] = U.hh2_1pt_2ft ;
            hp2.V_1ft_2np[1,dHH,k,l,.,.] = U.hh2_1ft_2np ;
            hp2.V_1ft_2ue[1,dHH,k,l,.,.] = U.hh2_1ft_2ue ;
            hp2.V_1ft_2pt[1,dHH,k,l,.,.] = U.hh2_1ft_2pt ;
            hp2.V_1ft_2ft[1,dHH,k,l,.,.] = U.hh2_1ft_2ft ;
 
            hp2.V_1np_2np_nk[1,dHH,k,l,.,.] = U.hh2_1np_2np_nk;
            hp2.V_1np_2ue_nk[1,dHH,k,l,.,.] = U.hh2_1np_2ue_nk;
            hp2.V_1np_2pt_nk[1,dHH,k,l,.,.] = U.hh2_1np_2pt_nk;
            hp2.V_1np_2ft_nk[1,dHH,k,l,.,.] = U.hh2_1np_2ft_nk;
            hp2.V_1ue_2np_nk[1,dHH,k,l,.,.] = U.hh2_1ue_2np_nk;
            hp2.V_1ue_2ue_nk[1,dHH,k,l,.,.] = U.hh2_1ue_2ue_nk;
            hp2.V_1ue_2pt_nk[1,dHH,k,l,.,.] = U.hh2_1ue_2pt_nk;
            hp2.V_1ue_2ft_nk[1,dHH,k,l,.,.] = U.hh2_1ue_2ft_nk;
            hp2.V_1pt_2np_nk[1,dHH,k,l,.,.] = U.hh2_1pt_2np_nk;
            hp2.V_1pt_2ue_nk[1,dHH,k,l,.,.] = U.hh2_1pt_2ue_nk;
            hp2.V_1pt_2pt_nk[1,dHH,k,l,.,.] = U.hh2_1pt_2pt_nk;
            hp2.V_1pt_2ft_nk[1,dHH,k,l,.,.] = U.hh2_1pt_2ft_nk;
            hp2.V_1ft_2np_nk[1,dHH,k,l,.,.] = U.hh2_1ft_2np_nk;
            hp2.V_1ft_2ue_nk[1,dHH,k,l,.,.] = U.hh2_1ft_2ue_nk;
            hp2.V_1ft_2pt_nk[1,dHH,k,l,.,.] = U.hh2_1ft_2pt_nk;
            hp2.V_1ft_2ft_nk[1,dHH,k,l,.,.] = U.hh2_1ft_2ft_nk;
 
            fixed_point1(V,U,PDF_ln);
 
            //Value function
            V_nl.V_1np_2np[1,dHH,k,l,.,.] = New_V.V_1np_2np ;
            V_nl.V_1np_2ue[1,dHH,k,l,.,.] = New_V.V_1np_2ue ;
            V_nl.V_1np_2pt[1,dHH,k,l,.,.] = New_V.V_1np_2pt ;
            V_nl.V_1np_2ft[1,dHH,k,l,.,.] = New_V.V_1np_2ft ;
            V_nl.V_1ue_2np[1,dHH,k,l,.,.] = New_V.V_1ue_2np ;
            V_nl.V_1ue_2ue[1,dHH,k,l,.,.] = New_V.V_1ue_2ue ;
            V_nl.V_1ue_2pt[1,dHH,k,l,.,.] = New_V.V_1ue_2pt ;
            V_nl.V_1ue_2ft[1,dHH,k,l,.,.] = New_V.V_1ue_2ft ;
            V_nl.V_1pt_2np[1,dHH,k,l,.,.] = New_V.V_1pt_2np ;
            V_nl.V_1pt_2ue[1,dHH,k,l,.,.] = New_V.V_1pt_2ue ;
            V_nl.V_1pt_2pt[1,dHH,k,l,.,.] = New_V.V_1pt_2pt ;
            V_nl.V_1pt_2ft[1,dHH,k,l,.,.] = New_V.V_1pt_2ft ;
            V_nl.V_1ft_2np[1,dHH,k,l,.,.] = New_V.V_1ft_2np ;
            V_nl.V_1ft_2ue[1,dHH,k,l,.,.] = New_V.V_1ft_2ue ;
            V_nl.V_1ft_2pt[1,dHH,k,l,.,.] = New_V.V_1ft_2pt ;
            V_nl.V_1ft_2ft[1,dHH,k,l,.,.] = New_V.V_1ft_2ft ;
 
            V_nl.V_1np_2np_nk[1,dHH,k,l,.,.] = New_V.V_1np_2np_nk;
            V_nl.V_1np_2ue_nk[1,dHH,k,l,.,.] = New_V.V_1np_2ue_nk;
            V_nl.V_1np_2pt_nk[1,dHH,k,l,.,.] = New_V.V_1np_2pt_nk;
            V_nl.V_1np_2ft_nk[1,dHH,k,l,.,.] = New_V.V_1np_2ft_nk;
            V_nl.V_1ue_2np_nk[1,dHH,k,l,.,.] = New_V.V_1ue_2np_nk;
            V_nl.V_1ue_2ue_nk[1,dHH,k,l,.,.] = New_V.V_1ue_2ue_nk;
            V_nl.V_1ue_2pt_nk[1,dHH,k,l,.,.] = New_V.V_1ue_2pt_nk;
            V_nl.V_1ue_2ft_nk[1,dHH,k,l,.,.] = New_V.V_1ue_2ft_nk;
            V_nl.V_1pt_2np_nk[1,dHH,k,l,.,.] = New_V.V_1pt_2np_nk;
            V_nl.V_1pt_2ue_nk[1,dHH,k,l,.,.] = New_V.V_1pt_2ue_nk;
            V_nl.V_1pt_2pt_nk[1,dHH,k,l,.,.] = New_V.V_1pt_2pt_nk;
            V_nl.V_1pt_2ft_nk[1,dHH,k,l,.,.] = New_V.V_1pt_2ft_nk;
            V_nl.V_1ft_2np_nk[1,dHH,k,l,.,.] = New_V.V_1ft_2np_nk;
            V_nl.V_1ft_2ue_nk[1,dHH,k,l,.,.] = New_V.V_1ft_2ue_nk;
            V_nl.V_1ft_2pt_nk[1,dHH,k,l,.,.] = New_V.V_1ft_2pt_nk;
            V_nl.V_1ft_2ft_nk[1,dHH,k,l,.,.] = New_V.V_1ft_2ft_nk;
 
            time_out=hsec;
            //print "NL-1 Fixed Point Time" (time_out-time_in)/100 "Iter" ii_fxdpt1 dHH hh_om;
 
            // Non-Labor Income 2
            time_in=hsec;
 
            V.V_1np_2np = arraytomat(V_nl.V_1np_2np[2,dHH,k,l,.,.]) ;
            V.V_1ue_2np = arraytomat(V_nl.V_1ue_2np[2,dHH,k,l,.,.]) ;
            V.V_1np_2ue = arraytomat(V_nl.V_1np_2ue[2,dHH,k,l,.,.]) ;
            V.V_1ue_2ue = arraytomat(V_nl.V_1ue_2ue[2,dHH,k,l,.,.]) ;
            V.V_1np_2pt = arraytomat(V_nl.V_1np_2pt[2,dHH,k,l,.,.]) ;
            V.V_1ue_2pt = arraytomat(V_nl.V_1ue_2pt[2,dHH,k,l,.,.]) ;
            V.V_1np_2ft = arraytomat(V_nl.V_1np_2ft[2,dHH,k,l,.,.]) ;
            V.V_1ue_2ft = arraytomat(V_nl.V_1ue_2ft[2,dHH,k,l,.,.]) ;
 
            V.V_1pt_2np = arraytomat(V_nl.V_1pt_2np[2,dHH,k,l,.,.]) ;
            V.V_1ft_2np = arraytomat(V_nl.V_1ft_2np[2,dHH,k,l,.,.]) ;
            V.V_1pt_2ue = arraytomat(V_nl.V_1pt_2ue[2,dHH,k,l,.,.]) ;
            V.V_1ft_2ue = arraytomat(V_nl.V_1ft_2ue[2,dHH,k,l,.,.]) ;
            V.V_1pt_2pt = arraytomat(V_nl.V_1pt_2pt[2,dHH,k,l,.,.]) ;
            V.V_1ft_2pt = arraytomat(V_nl.V_1ft_2pt[2,dHH,k,l,.,.]) ;
            V.V_1pt_2ft = arraytomat(V_nl.V_1pt_2ft[2,dHH,k,l,.,.]) ;
            V.V_1ft_2ft = arraytomat(V_nl.V_1ft_2ft[2,dHH,k,l,.,.]) ;
 
            V.V_1np_2np_nk = arraytomat(V_nl.V_1np_2np_nk[2,dHH,k,l,.,.]);
            V.V_1ue_2np_nk = arraytomat(V_nl.V_1ue_2np_nk[2,dHH,k,l,.,.]);
            V.V_1np_2ue_nk = arraytomat(V_nl.V_1np_2ue_nk[2,dHH,k,l,.,.]);
            V.V_1ue_2ue_nk = arraytomat(V_nl.V_1ue_2ue_nk[2,dHH,k,l,.,.]);
            V.V_1np_2pt_nk = arraytomat(V_nl.V_1np_2pt_nk[2,dHH,k,l,.,.]);
            V.V_1ue_2pt_nk = arraytomat(V_nl.V_1ue_2pt_nk[2,dHH,k,l,.,.]);
            V.V_1np_2ft_nk = arraytomat(V_nl.V_1np_2ft_nk[2,dHH,k,l,.,.]);
            V.V_1ue_2ft_nk = arraytomat(V_nl.V_1ue_2ft_nk[2,dHH,k,l,.,.]);
 
            V.V_1pt_2np_nk = arraytomat(V_nl.V_1pt_2np_nk[2,dHH,k,l,.,.]);
            V.V_1ft_2np_nk = arraytomat(V_nl.V_1ft_2np_nk[2,dHH,k,l,.,.]);
            V.V_1pt_2ue_nk = arraytomat(V_nl.V_1pt_2ue_nk[2,dHH,k,l,.,.]);
            V.V_1ft_2ue_nk = arraytomat(V_nl.V_1ft_2ue_nk[2,dHH,k,l,.,.]);
            V.V_1pt_2pt_nk = arraytomat(V_nl.V_1pt_2pt_nk[2,dHH,k,l,.,.]);
            V.V_1ft_2pt_nk = arraytomat(V_nl.V_1ft_2pt_nk[2,dHH,k,l,.,.]);
            V.V_1pt_2ft_nk = arraytomat(V_nl.V_1pt_2ft_nk[2,dHH,k,l,.,.]);
            V.V_1ft_2ft_nk = arraytomat(V_nl.V_1ft_2ft_nk[2,dHH,k,l,.,.]);
 
            U = hh_utility(nonlabInc2,dHH);
 
            hp1.V_1np_2np[2,dHH,k,l,.,.] = U.hh1_1np_2np ;
            hp1.V_1np_2ue[2,dHH,k,l,.,.] = U.hh1_1np_2ue ;
            hp1.V_1np_2pt[2,dHH,k,l,.,.] = U.hh1_1np_2pt ;
            hp1.V_1np_2ft[2,dHH,k,l,.,.] = U.hh1_1np_2ft ;
            hp1.V_1ue_2np[2,dHH,k,l,.,.] = U.hh1_1ue_2np ;
            hp1.V_1ue_2ue[2,dHH,k,l,.,.] = U.hh1_1ue_2ue ;
            hp1.V_1ue_2pt[2,dHH,k,l,.,.] = U.hh1_1ue_2pt ;
            hp1.V_1ue_2ft[2,dHH,k,l,.,.] = U.hh1_1ue_2ft ;
            hp1.V_1pt_2np[2,dHH,k,l,.,.] = U.hh1_1pt_2np ;
            hp1.V_1pt_2ue[2,dHH,k,l,.,.] = U.hh1_1pt_2ue ;
            hp1.V_1pt_2pt[2,dHH,k,l,.,.] = U.hh1_1pt_2pt ;
            hp1.V_1pt_2ft[2,dHH,k,l,.,.] = U.hh1_1pt_2ft ;
            hp1.V_1ft_2np[2,dHH,k,l,.,.] = U.hh1_1ft_2np ;
            hp1.V_1ft_2ue[2,dHH,k,l,.,.] = U.hh1_1ft_2ue ;
            hp1.V_1ft_2pt[2,dHH,k,l,.,.] = U.hh1_1ft_2pt ;
            hp1.V_1ft_2ft[2,dHH,k,l,.,.] = U.hh1_1ft_2ft ;
 
            hp1.V_1np_2np_nk[2,dHH,k,l,.,.] = U.hh1_1np_2np_nk;
            hp1.V_1np_2ue_nk[2,dHH,k,l,.,.] = U.hh1_1np_2ue_nk;
            hp1.V_1np_2pt_nk[2,dHH,k,l,.,.] = U.hh1_1np_2pt_nk;
            hp1.V_1np_2ft_nk[2,dHH,k,l,.,.] = U.hh1_1np_2ft_nk;
            hp1.V_1ue_2np_nk[2,dHH,k,l,.,.] = U.hh1_1ue_2np_nk;
            hp1.V_1ue_2ue_nk[2,dHH,k,l,.,.] = U.hh1_1ue_2ue_nk;
            hp1.V_1ue_2pt_nk[2,dHH,k,l,.,.] = U.hh1_1ue_2pt_nk;
            hp1.V_1ue_2ft_nk[2,dHH,k,l,.,.] = U.hh1_1ue_2ft_nk;
            hp1.V_1pt_2np_nk[2,dHH,k,l,.,.] = U.hh1_1pt_2np_nk;
            hp1.V_1pt_2ue_nk[2,dHH,k,l,.,.] = U.hh1_1pt_2ue_nk;
            hp1.V_1pt_2pt_nk[2,dHH,k,l,.,.] = U.hh1_1pt_2pt_nk;
            hp1.V_1pt_2ft_nk[2,dHH,k,l,.,.] = U.hh1_1pt_2ft_nk;
            hp1.V_1ft_2np_nk[2,dHH,k,l,.,.] = U.hh1_1ft_2np_nk;
            hp1.V_1ft_2ue_nk[2,dHH,k,l,.,.] = U.hh1_1ft_2ue_nk;
            hp1.V_1ft_2pt_nk[2,dHH,k,l,.,.] = U.hh1_1ft_2pt_nk;
            hp1.V_1ft_2ft_nk[2,dHH,k,l,.,.] = U.hh1_1ft_2ft_nk;
 
            hp2.V_1np_2np[2,dHH,k,l,.,.] = U.hh2_1np_2np ;
            hp2.V_1np_2ue[2,dHH,k,l,.,.] = U.hh2_1np_2ue ;
            hp2.V_1np_2pt[2,dHH,k,l,.,.] = U.hh2_1np_2pt ;
            hp2.V_1np_2ft[2,dHH,k,l,.,.] = U.hh2_1np_2ft ;
            hp2.V_1ue_2np[2,dHH,k,l,.,.] = U.hh2_1ue_2np ;
            hp2.V_1ue_2ue[2,dHH,k,l,.,.] = U.hh2_1ue_2ue ;
            hp2.V_1ue_2pt[2,dHH,k,l,.,.] = U.hh2_1ue_2pt ;
            hp2.V_1ue_2ft[2,dHH,k,l,.,.] = U.hh2_1ue_2ft ;
            hp2.V_1pt_2np[2,dHH,k,l,.,.] = U.hh2_1pt_2np ;
            hp2.V_1pt_2ue[2,dHH,k,l,.,.] = U.hh2_1pt_2ue ;
            hp2.V_1pt_2pt[2,dHH,k,l,.,.] = U.hh2_1pt_2pt ;
            hp2.V_1pt_2ft[2,dHH,k,l,.,.] = U.hh2_1pt_2ft ;
            hp2.V_1ft_2np[2,dHH,k,l,.,.] = U.hh2_1ft_2np ;
            hp2.V_1ft_2ue[2,dHH,k,l,.,.] = U.hh2_1ft_2ue ;
            hp2.V_1ft_2pt[2,dHH,k,l,.,.] = U.hh2_1ft_2pt ;
            hp2.V_1ft_2ft[2,dHH,k,l,.,.] = U.hh2_1ft_2ft ;
 
            hp2.V_1np_2np_nk[2,dHH,k,l,.,.] = U.hh2_1np_2np_nk;
            hp2.V_1np_2ue_nk[2,dHH,k,l,.,.] = U.hh2_1np_2ue_nk;
            hp2.V_1np_2pt_nk[2,dHH,k,l,.,.] = U.hh2_1np_2pt_nk;
            hp2.V_1np_2ft_nk[2,dHH,k,l,.,.] = U.hh2_1np_2ft_nk;
            hp2.V_1ue_2np_nk[2,dHH,k,l,.,.] = U.hh2_1ue_2np_nk;
            hp2.V_1ue_2ue_nk[2,dHH,k,l,.,.] = U.hh2_1ue_2ue_nk;
            hp2.V_1ue_2pt_nk[2,dHH,k,l,.,.] = U.hh2_1ue_2pt_nk;
            hp2.V_1ue_2ft_nk[2,dHH,k,l,.,.] = U.hh2_1ue_2ft_nk;
            hp2.V_1pt_2np_nk[2,dHH,k,l,.,.] = U.hh2_1pt_2np_nk;
            hp2.V_1pt_2ue_nk[2,dHH,k,l,.,.] = U.hh2_1pt_2ue_nk;
            hp2.V_1pt_2pt_nk[2,dHH,k,l,.,.] = U.hh2_1pt_2pt_nk;
            hp2.V_1pt_2ft_nk[2,dHH,k,l,.,.] = U.hh2_1pt_2ft_nk;
            hp2.V_1ft_2np_nk[2,dHH,k,l,.,.] = U.hh2_1ft_2np_nk;
            hp2.V_1ft_2ue_nk[2,dHH,k,l,.,.] = U.hh2_1ft_2ue_nk;
            hp2.V_1ft_2pt_nk[2,dHH,k,l,.,.] = U.hh2_1ft_2pt_nk;
            hp2.V_1ft_2ft_nk[2,dHH,k,l,.,.] = U.hh2_1ft_2ft_nk;
 
            fixed_point1(V,U,PDF_ln);
 
            V_nl.V_1np_2np[2,dHH,k,l,.,.] = New_V.V_1np_2np ;
            V_nl.V_1np_2ue[2,dHH,k,l,.,.] = New_V.V_1np_2ue ;
            V_nl.V_1np_2pt[2,dHH,k,l,.,.] = New_V.V_1np_2pt ;
            V_nl.V_1np_2ft[2,dHH,k,l,.,.] = New_V.V_1np_2ft ;
            V_nl.V_1ue_2np[2,dHH,k,l,.,.] = New_V.V_1ue_2np ;
            V_nl.V_1ue_2ue[2,dHH,k,l,.,.] = New_V.V_1ue_2ue ;
            V_nl.V_1ue_2pt[2,dHH,k,l,.,.] = New_V.V_1ue_2pt ;
            V_nl.V_1ue_2ft[2,dHH,k,l,.,.] = New_V.V_1ue_2ft ;
            V_nl.V_1pt_2np[2,dHH,k,l,.,.] = New_V.V_1pt_2np ;
            V_nl.V_1pt_2ue[2,dHH,k,l,.,.] = New_V.V_1pt_2ue ;
            V_nl.V_1pt_2pt[2,dHH,k,l,.,.] = New_V.V_1pt_2pt ;
            V_nl.V_1pt_2ft[2,dHH,k,l,.,.] = New_V.V_1pt_2ft ;
            V_nl.V_1ft_2np[2,dHH,k,l,.,.] = New_V.V_1ft_2np ;
            V_nl.V_1ft_2ue[2,dHH,k,l,.,.] = New_V.V_1ft_2ue ;
            V_nl.V_1ft_2pt[2,dHH,k,l,.,.] = New_V.V_1ft_2pt ;
            V_nl.V_1ft_2ft[2,dHH,k,l,.,.] = New_V.V_1ft_2ft ;
 
            V_nl.V_1np_2np_nk[2,dHH,k,l,.,.] = New_V.V_1np_2np_nk;
            V_nl.V_1np_2ue_nk[2,dHH,k,l,.,.] = New_V.V_1np_2ue_nk;
            V_nl.V_1np_2pt_nk[2,dHH,k,l,.,.] = New_V.V_1np_2pt_nk;
            V_nl.V_1np_2ft_nk[2,dHH,k,l,.,.] = New_V.V_1np_2ft_nk;
            V_nl.V_1ue_2np_nk[2,dHH,k,l,.,.] = New_V.V_1ue_2np_nk;
            V_nl.V_1ue_2ue_nk[2,dHH,k,l,.,.] = New_V.V_1ue_2ue_nk;
            V_nl.V_1ue_2pt_nk[2,dHH,k,l,.,.] = New_V.V_1ue_2pt_nk;
            V_nl.V_1ue_2ft_nk[2,dHH,k,l,.,.] = New_V.V_1ue_2ft_nk;
            V_nl.V_1pt_2np_nk[2,dHH,k,l,.,.] = New_V.V_1pt_2np_nk;
            V_nl.V_1pt_2ue_nk[2,dHH,k,l,.,.] = New_V.V_1pt_2ue_nk;
            V_nl.V_1pt_2pt_nk[2,dHH,k,l,.,.] = New_V.V_1pt_2pt_nk;
            V_nl.V_1pt_2ft_nk[2,dHH,k,l,.,.] = New_V.V_1pt_2ft_nk;
            V_nl.V_1ft_2np_nk[2,dHH,k,l,.,.] = New_V.V_1ft_2np_nk;
            V_nl.V_1ft_2ue_nk[2,dHH,k,l,.,.] = New_V.V_1ft_2ue_nk;
            V_nl.V_1ft_2pt_nk[2,dHH,k,l,.,.] = New_V.V_1ft_2pt_nk;
            V_nl.V_1ft_2ft_nk[2,dHH,k,l,.,.] = New_V.V_1ft_2ft_nk;
 
            time_out = hsec ;
            time_out_s = hsec ;
            //print "NL-2 Fixed Point Time" (time_out-time_in)/100 "Iter" ii_fxdpt1 dHH hh_om;
 
            //print "Total Time per Search Cost:" (time_out_s-time_in_s)/100;
        endfor;
    endfor;
endfor;
/* Initialization of Variables */
spouses_array=arrayinit(N_Agents | column_matrices_max | 15,0);
duration_matrix = zeros(N_Agents,column_matrices_max) ;
//order state wage_m wage_f hp_m hp_m kids value 8 9 10 11 duration
 
ii_total = ones(N_agents,1).*2 ; /* spell_counter *///2=labor shocks needs to be updated when kids shock arrived 1st (compared with t-1)
ii_total_k = ones(N_agents,1) ; /* kids spell_counter */
ii_labor_m = ones(N_agents,1) ; /* firing shock males*/
ii_labor_f = ones(N_agents,1) ; /* firing shock females*/
 
/* Simulation */
//Threadfor i (1,N_agents,1) ;//i=1;
for i (1,N_agents,1) ;//i=1;
 
    ii_wage_m = 1;
    ii_wage_f = 1;
    ii_total_i = 2;//always start at 2
    ii_total_ki = 1;
    ii_next_state = 1;
    ind_kids = 0; /*indicator if the kids shock had happened*/
    next_state = 1; /*formerly zeros(N_agents,column_matrices_max)*/
    labor_time = 0;
    n_sta = 0;
    state = 0;
    wage_m_new = 1;
    wage_f_new = 1;
    time_m_emp = 1;
    time_m_exo = 1;
    time_m_j2j = 1;
    time_f_emp = 1;
    time_f_exo = 1;
    time_f_j2j = 1;
    NewVF=0;
 
    spouses_matrix =zeros(column_matrices_max,15) ;
 
    //Random matrices of shocks per individual
    zeta_nkids2=zeta_nkids[i,.];
    zeta_ykids2=zeta_ykids[i,.];
 
    zeta_eta_m2=zeta_eta_m[i,.];
    zeta_j2j_m2=zeta_j2j_m[i,.];
    zeta_u2e_m2=zeta_u2e_m[i,.];
    zeta_eta_f2=zeta_eta_f[i,.];
    zeta_j2j_f2=zeta_j2j_f[i,.];
    zeta_u2e_f2=zeta_u2e_f[i,.];
 
    zeta_hours_m2=zeta_hours_m[i,.];
    zeta_hours_f2=zeta_hours_f[i,.];
    zeta_wage_m2=zeta_wage_m[i,.];
    zeta_wage_f2=zeta_wage_f[i,.];
 
    // now see how long it takes for each spouse to get exogenously dismissed from current job
    time_m_exo=-ln(1-zeta_eta_m2[ii_wage_m ])/eta1_pt;
    // now see how long it takes each spouse alternative job offer, whether pt or ft
    time_m_j2j=-ln(1-zeta_j2j_m2[ii_wage_m ])/lam_e1;
    // now see how long it takes for each spouse to get exogenously dismissed from current job
    time_f_exo=-ln(1-zeta_eta_f2[ii_wage_f ])/eta2_pt;
    // now see how long it takes each spouse alternative job offer, whether pt or ft
    time_f_j2j=-ln(1-zeta_j2j_f2[ii_wage_f ])/lam_e2;
 
    //time_m_emp= -ln(1-zeta_u2e_m2[ii_wage_m])/lam_n1;
    //time_f_emp= -ln(1-zeta_u2e_f2[ii_wage_f])/lam_n2;
 
    //Initial conditions
    //Labor market initial state: pt-pt:wages==3
    next_state=6 ;
    wage_m = 3;
    wage_f = 3;
    w_ind_m=minindc( miss(maxc( (w1_vector~wage_m*ones(N_w1 ,1) )') , wage_m) ) ;
    w_ind_f=minindc( miss(maxc( (w2_vector~wage_f*ones(N_w1 ,1) )') , wage_f) ) ;
 
    // Search Cost Indexes here productivity types (k,l)
    k = simul_s1[i] ; //male
    l = simul_s2[i] ; //female
    // Children, Non Labor Income and Other members initial status
    ninos=kids_matrix1[i];
    simul_Nlab_i=simul_nonlabor_vec[i];
    d_HH_i =d13_matrix[i];
 
    if d_HH_i==1; // 1=Other members 2=No Other members
        gam_nk2=gam_nkd;
        gam_yk2=gam_ykd;
    else;
        gam_nk2=gam_nk ;
        gam_yk2=gam_yk ;
    endif;
 
    if ninos==1; // 1=Have children 2=Not have children
        zeta_kids1=zeta_nkids2;
        gam_kids=gam_nk2;
    else;
        zeta_kids1=zeta_ykids2;
        gam_kids=gam_yk2;
    endif;
 
    time_kids=-ln(1-zeta_kids1[ii_total_ki])/ gam_kids ;
    labor_time=time_m_exo|time_m_j2j|time_f_exo|time_f_j2j ;
 
    mat_vec(ninos,d_HH_i);
    //threadEndFor;
    VF_value=VF_1pt_2pt[w_ind_m,w_ind_f ];
    //print k l VF_1np_2np;
    // This indicator tell us if a kids shock was the last one
    ind_kids=0;
    // This indicator tell us if a labor shock was the last one
    switch_shock_m=0;
    switch_shock_f=0;
 
    do until sumc(spouses_matrix[.,12]).>=length_labor_career ;
 
        //wage_m wage_f VF_value state_matrix[i,ii_total_i] HHprod_matrix[i,ii_total_i]
        //ii_total_i ii_wage_m ii_wage_f ii_labor_m[i] ii_labor_f[i];
        //if i<10;
        //print "here0" i spouses_matrix[ii_total_i-1,12] next_state next_hh_prod VF_value wage_m wage_f ii_next_state NewVF wage_m_new wage_f_new ind_kids ninos gam_kids ii_total_i ii_wage_m ii_wage_f ii_total_ki switch_shock_m switch_shock_f zeta_hours_m2[ii_wage_m ].>prob1_pt time_kids labor_time' "a" time_m_emp time_m_exo time_m_j2j time_f_emp time_f_exo time_f_j2j;
        //endif;
        /* States:
        1 -> Husband UE - Wife UE 1
        2 -> Husband PT - Wife UE 2
        3 -> Husband FT - Wife UE 3
        4 -> Husband UE - Wife PT 4
        5 -> Husband UE - Wife FT 5
        6 -> Husband PT - Wife PT 6
        7 -> Husband FT - Wife PT 7
        8 -> Husband PT - Wife FT 8
        9 -> Husband FT - Wife FT 9
        10 -> Husband PT - Wife NP 10
        11 -> Husband FT - Wife NP 11
        12 -> Husband NP - Wife PT 12
        13 -> Husband NP - Wife FT 13
        14 -> Husband UE - Wife NP 14
        15 -> Husband NP - Wife UE 15
        16 -> Husband NP - Wife NP 16
        */
 
        /********* STATE 1 Husband UE - Wife UE***************/
        if next_state.==1 ;
            spouses_matrix[ii_total_i,1:7]=next_state~wage_m~wage_f~hp1_1ue_2ue~hp2_1ue_2ue~ninos~VF_value;
 
            if switch_shock_m==1; //not resample males shock
                /* now see how long it takes each spouse to get a job offer, whether pt or ft */
                time_m_emp=time_m_emp -switch_shock_m*spouses_matrix[ii_total_i-1,12];
            else;
                /* now see how long it takes each spouse to get a job offer, whether pt or ft */
                time_m_emp= -ln(1-zeta_u2e_m2[ii_wage_m])/lam_n1;
            endif;
            if switch_shock_f==1; //not resample females shock
                /* now see how long it takes each spouse to get a job offer, whether pt or ft */
                time_f_emp= time_f_emp -switch_shock_f*spouses_matrix[ii_total_i-1,12];
            else;
                /* now see how long it takes each spouse to get a job offer, whether pt or ft */
                time_f_emp= -ln(1-zeta_u2e_f2[ii_wage_f])/lam_n2;
            endif;
 
            //endif;
            if ind_kids.==0; //no kids shock in las period
                time_kids=time_kids-spouses_matrix[ii_total_i-1,12];
            endif;
            /*Labor time*/
            labor_time=time_m_emp|time_f_emp ;
 
            if minc(labor_time)<time_kids;
                spouses_matrix[ii_total_i,12]=minc(labor_time);
 
                if minindc(labor_time)==1; //time_m_emp|time_f_emp ;
                    /*Male job offer*/
                    if zeta_hours_m2[ii_wage_m ].<=prob1_pt ; //part time offer
                        pf_ft_m=1; //pt=1 ft=2
                        w_indPot_m=minindc( miss(maxc( (P_M_pt_cum~zeta_wage_m2[ii_wage_m]*ones(N_w1 ,1) )') , zeta_wage_m2[ii_wage_m]) ) ;
                        wage_m_new = w1_vector[w_indPot_m] ;
                        VF_dec=
                            VF_1pt_2ue[w_indPot_m] ~
                            VF_1pt_2np[w_indPot_m] ;
 
                        ind_dec=maxindc(VF_dec');
                        NewVF= VF_dec[ind_dec];
                        string VFpot_str = {"1pt_2ue_" "1pt_2np_"};
 
                    elseif zeta_hours_m2[ii_wage_m ].>prob1_pt ; //full time offer
                        pf_ft_m=2; //pt=1 ft=2
                        w_indPot_m=minindc( miss(maxc( (P_M_ft_cum~zeta_wage_m2[ii_wage_m]*ones(N_w1 ,1) )') , zeta_wage_m2[ii_wage_m]) ) ;
                        wage_m_new = w1_vector[w_indPot_m] ;
                        VF_dec=
                            VF_1ft_2ue[w_indPot_m] ~
                            VF_1ft_2np[w_indPot_m] ;
 
                        ind_dec=maxindc(VF_dec');
                        NewVF= VF_dec[ind_dec];
                        string VFpot_str = {"1ft_2ue_" "1ft_2np_"};
 
                    endif;
                    /*To choose the labor states decision line*/
                    ii_next_state=ind_dec;
                    //new shock counter
                    ii_total_i=ii_total_i+1 ;
                    //new wage shock
                    ii_wage_m=ii_wage_m+1;
 
                    switch_shock_m=0; //It was a Male Shock
                    switch_shock_f=1; //Female initialy not resample (Only if endogenous quit)
 
                    if NewVF<=VF_value; //remains in the same state
                        VF_value = VF_value ;
                        wage_m = 0 ;
                        wage_f = 0 ;
 
                    elseif NewVF>VF_value; //Switch state
                        VF_value = NewVF ;
 
                        if zeta_hours_m2[ii_wage_m-1 ].<=prob1_pt ; //part time offer //"1pt_2ue_" "1pt_2pt_" "1pt_2np_"
                            if ii_next_state==1; //PT_UE
                                next_state=2;
                                wage_m = wage_m_new ;
                                wage_f = 0 ;
                                w_ind_m=w_indPot_m;
 
                            elseif ii_next_state==2; //PT_NP
                                next_state=10;
                                wage_m = wage_m_new ;
                                wage_f = 0 ;
                                w_ind_m=w_indPot_m;
                                switch_shock_f=0; //Female endogenous quit
                            endif;
 
                        elseif zeta_hours_m2[ii_wage_m-1 ].>prob1_pt ; //full time offer //"1ft_2ue_" "1ft_2pt_" "1ft_2np_""
                            if ii_next_state==1; //FT_UE
                                next_state=3;
                                wage_m = wage_m_new ;
                                wage_f = 0 ;
                                w_ind_m=w_indPot_m;
 
                            elseif ii_next_state==2; //FT_NP
                                next_state=11;
                                wage_m = wage_m_new ;
                                wage_f = 0 ;
                                w_ind_m=w_indPot_m;
                                switch_shock_f=0; //Female endogenous quit
                            endif;
                        endif;
                    endif;
 
                elseif minindc(labor_time)==2; //time_m_emp|time_f_emp ;
                    /*Female job offer*/
                    if zeta_hours_f2[ii_wage_f ].<=prob2_pt ; //part time offer
                        pf_ft_f=1; //pt=1 ft=2
                        w_indPot_f=minindc( miss(maxc( (P_W_pt_cum~zeta_wage_f2[ii_wage_f]*ones(N_w1 ,1) )') , zeta_wage_f2[ii_wage_f]) ) ;
                        wage_f_new = w2_vector[w_indPot_f] ;
                        VF_dec=
                            VF_1ue_2pt[w_indPot_f] ~
                            VF_1np_2pt[w_indPot_f] ;
 
                        ind_dec=maxindc(VF_dec');
                        NewVF= VF_dec[ind_dec];
                        string VFpot_str = {"1ue_2pt_" "1np_2pt_"};
 
                    elseif zeta_hours_f2[ii_wage_f ].>prob2_pt ; //full time offer
                        pf_ft_f=2; //pt=1 ft=2
                        w_indPot_f=minindc( miss(maxc( (P_W_ft_cum~zeta_wage_f2[ii_wage_f]*ones(N_w1 ,1) )') , zeta_wage_f2[ii_wage_f]) ) ;
                        wage_f_new = w2_vector[w_indPot_f] ;
                        VF_dec=
                            VF_1ue_2ft[w_indPot_f] ~
                            VF_1np_2ft[w_indPot_f] ;
 
                        ind_dec=maxindc(VF_dec');
                        NewVF= VF_dec[ind_dec];
                        string VFpot_str = {"1ue_2ft_" "1np_2ft_"};
 
                    endif;
 
                    /*To choose the labor states decision line*/
                    ii_next_state=ind_dec;
                    //new shock counter
                    ii_total_i=ii_total_i+1 ;
                    //new wage shock
                    ii_wage_f=ii_wage_f+1;
 
                    switch_shock_m=1; //Male initialy not resample (Only if endogenous quit)
                    switch_shock_f=0; //It was a Female Shock
 
                    if NewVF<=VF_value; //remains in the same state
                        VF_value = VF_value ;
                        wage_m = 0 ;
                        wage_f = 0 ;
 
                    elseif NewVF>VF_value; //Switch state
                        VF_value = NewVF ;
 
                        if zeta_hours_f2[ii_wage_f-1 ].<=prob2_pt ; //part time offer //"1ue_2pt_" "1pt_2pt_" "1np_2pt_"
                            if ii_next_state==1; //UE_PT
                                next_state=4;
                                wage_m = 0 ;
                                wage_f = wage_f_new ;
                                w_ind_f=w_indPot_f;
 
                            elseif ii_next_state==2; //NP_PT
                                next_state=12;
                                wage_m = 0 ;
                                wage_f = wage_f_new ;
                                w_ind_f=w_indPot_f;
                                switch_shock_m=0; //Male endogenous quit
                            endif;
 
                        elseif zeta_hours_f2[ii_wage_f-1 ].>prob2_pt ; //full time offer //"1ue_2ft_" "1pt_2ft_" "1np_2ft_"
                            if ii_next_state==1; //UE_FT
                                next_state=5;
                                wage_m = 0 ;
                                wage_f = wage_f_new ;
                                w_ind_f=w_indPot_f;
 
                            elseif ii_next_state==2; //NP_FT
                                next_state=13;
                                wage_m = 0 ;
                                wage_f = wage_f_new ;
                                w_ind_f=w_indPot_f;
                                switch_shock_m=0; //Male endogenous quit
                            endif;
 
                        endif;
                    endif;
 
                endif;
                ind_kids=0;
 
            else;
                spouses_matrix[ii_total_i,12]=time_kids;
 
                if ninos==1; //if in the present the couple has children
                    ninos=2; //the couple lose children
                    mat_vec(ninos,d_HH_i); //Change of the value functions
                    zeta_kids1=zeta_ykids2; //the next shock will be have children again
                    gam_kids=gam_yk2; //the next gamma will be have children again
                else; //if in the present the couple has no children
                    ninos=1; //the couple gain children
                    mat_vec(ninos,d_HH_i); //Change of the value functions
                    zeta_kids1=zeta_nkids2; //the next shock will be have no children again
                    gam_kids=gam_nk2; //the next gamma will be have no children again
                endif;
 
                VF_dec=
                    VF_1ue_2ue ~
                    VF_1ue_2np ~
                    VF_1np_2ue ~
                    VF_1np_2np ;
 
                ind_dec=maxindc(VF_dec');
                NewVF= VF_dec[ind_dec];
 
                string VFpot_str = {"1ue_2ue_" "1ue_2np_" "1np_2ue_" "1np_2np_"};
 
                /*To choose the labor states decision line*/
                ii_next_state=ind_dec;
                //new termination shock
                ii_total_i=ii_total_i+1 ;
                //new kids counter
                ii_total_ki=ii_total_ki+1;
 
                VF_value = NewVF ;
 
                //"1ue_2ue_" "1pt_2ue_" "1ue_2pt_" "1pt_2pt_" "1pt_2np_" "1np_2pt_" "1ue_2np_" "1np_2ue_" "1np_2np_"
                if ii_next_state==1; //UE_UE
                    next_state=1;
                    switch_shock_m=1;
                    switch_shock_f=1;
                elseif ii_next_state==2; //UE_NP
                    next_state=14;
                    switch_shock_m=1;
                    switch_shock_f=0;
                elseif ii_next_state==3; //NP_UE
                    next_state=15;
                    switch_shock_m=0;
                    switch_shock_f=1;
                elseif ii_next_state==4; //NP_NP
                    next_state=16;
                    switch_shock_m=0;
                    switch_shock_f=0;
                endif;
 
                wage_m = 0 ;
                wage_f = 0 ;
 
                //see how long it takes for couple to change the kids status
                time_kids=-ln(1-zeta_kids1[ii_total_ki])/ gam_kids ;
                ind_kids=1;
 
            endif;
 
        endif ;
        @ Simulation per HH is Finished -
        LIFE of and agent is simulated @
 
        //if i<10;
        //print "here1" i spouses_matrix[ii_total_i-1,12] next_state next_hh_prod VF_value wage_m wage_f ii_next_state NewVF wage_m_new wage_f_new ind_kids ninos gam_kids ii_total_i ii_wage_m ii_wage_f ii_total_ki switch_shock_m switch_shock_f zeta_hours_m2[ii_wage_m ].>prob1_pt time_kids labor_time' "a" time_m_emp time_m_exo time_m_j2j time_f_emp time_f_exo time_f_j2j;
        //endif;
    endo ; @ for sumc(duration_matrix[i,.]').>=length_labor_career do loop @
    //kids_matrix1[i]=ninos; // CHECK LF: this is the kids status in which the hh ends its life: why we need it?
    spouses_array[i,.,.]=spouses_matrix;
endfor ; @ corresponds to the for that starts the simulation for i (1,n,1) @
//threadEndFor;

And sometimes I am calling the Matvec that is outside.

I have looked over the code a bit. I think that you can remove the global assignments entirely. For example around line 682 of my version of the code (from earlier in this thread) looks like this:

if ninos==1; //if in the present the couple has children
    ninos=2; //the couple lose children
    mat_vec(ninos,d_HH_i); //Change of the value functions
    zeta_kids1=zeta_ykids2; //the next shock will be have children again
    gam_kids=gam_yk2; //the next gamma will be have children again
else; //if in the present the couple has no children
    ninos=1; //the couple gain children
    mat_vec(ninos,d_HH_i); //Change of the value functions
    zeta_kids1=zeta_nkids2; //the next shock will be have no children again
    gam_kids=gam_nk2; //the next gamma will be have no children again
endif;
 
VF_dec=
    VF_1ue_2ue ~
    VF_1ue_2np ~
    VF_1np_2ue ~
    VF_1np_2np ;
 
ind_dec=maxindc(VF_dec');

In this section of the code the only result of the mat_vec global variable assignments which are used are VF_1ue_2ue, VF_1ue_2np, VF_1np_2ue, VF_1np_2np. They are used to create ind_dec. If we create a procedure like this:

proc (1) = computeIndDec(ninos, dHH);
    local VF_dec;
    if ninos == 1;
        VF_dec = V_nl.V_1ue_2ue[simul_Nlab_i,dHH,k,l,.,.] ~
            V_nl.V_1ue_2np[simul_Nlab_i,dHH,k,l,.,.] ~
            V_nl.V_1np_2ue[simul_Nlab_i,dHH,k,l,.,.] ~
            V_nl.V_1np_2np[simul_Nlab_i,dHH,k,l,.,.];
    else;
        VF_dec = V_nl.V_1ue_2ue_nk[simul_Nlab_i,dHH,k,l,.,.] ~
            V_nl.V_1ue_2np_nk[simul_Nlab_i,dHH,k,l,.,.] ~
            V_nl.V_1np_2ue_nk[simul_Nlab_i,dHH,k,l,.,.] ~
            V_nl.V_1np_2np_nk[simul_Nlab_i,dHH,k,l,.,.];
    endif;
    retp(maxindc(VF_dec'));
endp;

and then change the first chunk of code to this:

if ninos==1; //if in the present the couple has children
    ninos=2; //the couple lose children
    zeta_kids1=zeta_ykids2; //the next shock will be have children again
    gam_kids=gam_yk2; //the next gamma will be have children again
else; //if in the present the couple has no children
    ninos=1; //the couple gain children
    zeta_kids1=zeta_nkids2; //the next shock will be have no children again
    gam_kids=gam_nk2; //the next gamma will be have no children again
endif;
 
ind_dec = computeIndDec(ninos, d_HH_i);

then we never need to make the global variable assignments to VF_1ue_2ue, VF_1ue_2np, VF_1np_2ue, VF_1np_2np.

If I did not miss anything and that works, I think it will make your code a bit cleaner as well. Let us know if that makes sense, or the solution will not work for some reason.

Hello, thanks for your answer. In fact Matvec is changing the values of VF_1ue_2ue, VF_1ue_2np, VF_1np_2ue, VF_1np_2np, but also of VF_1pt_2ue, VF_1pt_2np, VF_1ft_2ue, VF_1np_2ft, and so on that I use several times in the code. After the part that I sent, are 5 lines of decision making using these matrices.

elseif minindc(labor_time)==2; //time_m_emp|time_f_emp ;
    /*Female job offer*/
    if zeta_hours_f2[ii_wage_f ].<=prob2_pt ; //part time offer
        pf_ft_f=1; //pt=1 ft=2
        w_indPot_f=minindc( miss(maxc( (P_W_pt_cum~zeta_wage_f2[ii_wage_f]*ones(N_w1 ,1) )') , zeta_wage_f2[ii_wage_f]) ) ;
        wage_f_new = w2_vector[w_indPot_f] ;
        VF_dec=
            VF_1ue_2pt[w_indPot_f] ~
            VF_1np_2pt[w_indPot_f] ;
 
        ind_dec=maxindc(VF_dec');
        NewVF= VF_dec[ind_dec];
        string VFpot_str = {"1ue_2pt_" "1np_2pt_"};
 
    elseif zeta_hours_f2[ii_wage_f ].>prob2_pt ; //full time offer
        pf_ft_f=2; //pt=1 ft=2
        w_indPot_f=minindc( miss(maxc( (P_W_ft_cum~zeta_wage_f2[ii_wage_f]*ones(N_w1 ,1) )') , zeta_wage_f2[ii_wage_f]) ) ;
        wage_f_new = w2_vector[w_indPot_f] ;
        VF_dec=
            VF_1ue_2ft[w_indPot_f] ~
            VF_1np_2ft[w_indPot_f] ;
 
        ind_dec=maxindc(VF_dec');
        NewVF= VF_dec[ind_dec];
        string VFpot_str = {"1ue_2ft_" "1np_2ft_"};
 
    endif;

Think of that in terms of how your decision making changes when you have children or you don't have children. In time your reaction to labor market shocks changes. Below I sent you a part where I call

VF_1np_2np = V_nl.V_1np_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1np_2ue = V_nl.V_1np_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1np_2pt = V_nl.V_1np_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1np_2ft = V_nl.V_1np_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1ue_2np = V_nl.V_1ue_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1ue_2ue = V_nl.V_1ue_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1ue_2pt = V_nl.V_1ue_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1ue_2ft = V_nl.V_1ue_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1pt_2np = V_nl.V_1pt_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1pt_2ue = V_nl.V_1pt_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1pt_2pt = V_nl.V_1pt_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1pt_2ft = V_nl.V_1pt_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1ft_2np = V_nl.V_1ft_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1ft_2ue = V_nl.V_1ft_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
VF_1ft_2pt = V_nl.V_1ft_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; VF_1ft_2ft = V_nl.V_1ft_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
 
hp1_1np_2np = hp1.V_1np_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1np_2ue = hp1.V_1np_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1np_2pt = hp1.V_1np_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1np_2ft = hp1.V_1np_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1ue_2np = hp1.V_1ue_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1ue_2ue = hp1.V_1ue_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1ue_2pt = hp1.V_1ue_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1ue_2ft = hp1.V_1ue_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1pt_2np = hp1.V_1pt_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1pt_2ue = hp1.V_1pt_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1pt_2pt = hp1.V_1pt_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1pt_2ft = hp1.V_1pt_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1ft_2np = hp1.V_1ft_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1ft_2ue = hp1.V_1ft_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp1_1ft_2pt = hp1.V_1ft_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp1_1ft_2ft = hp1.V_1ft_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
 
hp2_1np_2np = hp2.V_1np_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1np_2ue = hp2.V_1np_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1np_2pt = hp2.V_1np_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1np_2ft = hp2.V_1np_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1ue_2np = hp2.V_1ue_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1ue_2ue = hp2.V_1ue_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1ue_2pt = hp2.V_1ue_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1ue_2ft = hp2.V_1ue_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1pt_2np = hp2.V_1pt_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1pt_2ue = hp2.V_1pt_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1pt_2pt = hp2.V_1pt_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1pt_2ft = hp2.V_1pt_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1ft_2np = hp2.V_1ft_2np_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1ft_2ue = hp2.V_1ft_2ue_nk[simul_Nlab_i,dHH,k,l,.,.];
hp2_1ft_2pt = hp2.V_1ft_2pt_nk[simul_Nlab_i,dHH,k,l,.,.]; hp2_1ft_2ft = hp2.V_1ft_2ft_nk[simul_Nlab_i,dHH,k,l,.,.];

And In the thread I call some locals that assing some values to those outputs, what do you think?