Looping with mata function

[Alan Wang]

Hi everyone, I happen to a problem when writing loop with mata function. Here I want to use grid search approach to solve for a set of functions and then store the results back to data. Hope someone could help me figure out how to make it work. Thank you in advance for reading my long thread...

The initial set up is:

set obs 20000
gen id=_n
gen e=.
gen phi=.
gen N_e=100 // number of values for e
gen N_phi=200 // number of values for phi
gen step_e=0.01 // step of grids for e
gen step_phi=1 // step of grids for phi
gen error_code=.
gen z1=.
gen z2=.
gen z3=.
gen b1=.
gen b2=.
gen b3=.

Then for any pair of values e and phi, I can use the following code to do what I want.

clear mata
mata:

temp_e=st_data(1, "e")
temp_phi=st_data(1, "phi")

st_local("e", strofreal(temp_e))
st_local("phi", strofreal(temp_phi))


void function solvez(real colvector z, real colvector y)
{

y[1] = ((20000+z[1])*(1-0.2)/(1+1/`e'))*(1-(20000/(20000+z[1]))^(1+1/`e'))-`phi'-z[1]*(1-0.25)
y[2] = ((9000+z[2])*(1-0.2)/(1+1/`e'))*(1-(9000/(9000+z[2]))^(1+1/`e'))-`phi'-z[2]*(1-0.25)
y[3] = ((35000+z[3])*(1-0.25)/(1+1/`e'))*(1-(35000/(35000+z[3]))^(1+1/`e'))-`phi'-z[3]*(1-0.3)

}
S = solvenl_init()
solvenl_init_evaluator(S, &solvez())
solvenl_init_type(S, "zero")
solvenl_init_technique(S, "newton")
solvenl_init_startingvals(S, J(3,1,10))
solvenl_init_conv_maxiter(S, 100)
solvenl_init_iter_log(S, "off")
_solvenl_solve(S)
ERR=solvenl_result_error_code(S)
M=solvenl_result_values(S)
st_matrix("M", M)
st_matrix("ERR", ERR)
end

matrix list M
matrix list ERR


replace z1=M[1,1] if id==`num'
replace z2=M[2,1] if id==`num'
replace z3=M[3,1] if id==`num'

replace b1=(e*0.05/(1-0.2)*20000-z1)/50 if id==`num'
replace b2=(e*0.05/(1-0.2)*9000-z2)/50 if id==`num'
replace b3=(e*0.05/(1-0.25)*35000-z3)/500 if id==`num'

replace error_code=ERR[1,1] if id==`num'


But I need to write a loop to do it for all 20000 pair of e and phi. Here is what I have tried my best to write.

capture program drop solveprogram
program define solveprogram
args num var_e var_phi

mata: mywork("`num'", "`var_e'", "`var_phi'")

replace z1=M[1,1] if id==`num'
replace z2=M[2,1] if id==`num'
replace z3=M[3,1] if id==`num'

replace b1=(e*0.05/(1-0.2)*20000-z1)/50 if id==`num'
replace b2=(e*0.05/(1-0.2)*9000-z2)/50 if id==`num'
replace b3=(e*0.05/(1-0.25)*35000-z3)/500 if id==`num'

replace error_code=ERR[1,1] if id==`num'
end

clear mata
mata:

void mywork(string scalar num_value, string scalar e_name, string scalar phi_name)
{

real scalar e_temp, phi_temp

e_temp = st_data(num_value, e_name)
phi_temp = st_data(num_value, phi_name)


S = solvenl_init()
solvenl_init_argument(S, 1, e_temp)
solvenl_init_argument(S, 2, phi_temp)
solvenl_init_evaluator(S, &solvez())
solvenl_init_type(S, "zero")
solvenl_init_technique(S, "newton")
solvenl_init_startingvals(S, J(3,1,10))
solvenl_init_conv_maxiter(S, 100)
solvenl_init_iter_log(S, "off")
_solvenl_solve(S)
ERR=solvenl_result_error_code(S)
M=solvenl_result_values(S)
st_matrix("M", M)
st_matrix("ERR", ERR)

}

void function solvez(real colvector z, real scalar e, real scalar phi, real colvector y)
{

y[1] = ((20000+z[1])*(1-0.2)/(1+1/e))*(1-(20000/(20000+z[1]))^(1+1/e))-phi-z[1]*(1-0.25)
y[2] = ((9000+z[2])*(1-0.2)/(1+1/e))*(1-(9000/(9000+z[2]))^(1+1/e))-phi-z[2]*(1-0.25)
y[3] = ((35000+z[3])*(1-0.25)/(1+1/e))*(1-(35000/(35000+z[3]))^(1+1/e))-phi-z[3]*(1-0.3)

}

end


solveprogram 1 e phi // try a case before writing the loop




And it reports error:
st_data(): 3253 nonreal found where real required
mywork(): - function returned error
<istmt>: - function returned error


I cannot figure out how to correct the problem since I don't have any nonreal values in the data. Could you help me make the code work?

[Christophe Kolodziejczyk]

the first argument of st_data() should be a number and not a string. You define scalar num_value as string and passes to st_data(), where it expects a number

[Alan Wang]

Thanks Christophe, could you tell me what should I revise then? I tried to rewrite void mywork(string scalar num_value) as void mywork(real scalar num_value), and then it returns
mywork(): 3253 <tmp>[1,1] found where real required
<istmt>: - function returned error

How could I deal with this then?

By the way, the initial set up should also include the following value endowment for e and phi.

foreach seq_e of numlist 1(1)100 {
foreach seq_phi of numlist 1(1)200 {
qui replace e=0.01*`seq_e' if id==(`seq_e'-1)*(N_phi)+`seq_phi'
qui replace phi=`seq_phi' if id==(`seq_e'-1)*(N_phi)+`seq_phi'

}
}

[Alan Wang]

My problem is solved. The main change is to rewrite void function solvez(real colvector z, real scalar e, real scalar phi, real colvector y) as void function solvez(real colvector z, real colvector y, real scalar e, real scalar phi).

[Nick Cox]

Cross-posted at http://stackoverflow.com/questions/3...lvenl-function

Please note our cross-posting policy, which is that you should tell us about it.