Extracting results from MATA

[Bernardo Coelho]

Hi everyone,
I'm using the Stata command blp to run a mixed logit. The problem is I need some results that the command generates but doesn't report.
More specifically, I need the vectors s_hat_t and S_hat_t generated in the ado file. Problem is, altough I'm quite used to use Stata, this is the first time i'm having to deal with MATA. I tried some simple things like saving a database after the vectors are generated, but wasn't successful.
One detail that is making this harder is that the vectors I want are generated within a void that uses a second command.
As my database has 5564 markets, the ideal would be that I could save this vector as a variable (the vector contains a value for each market, thus has the same size as my database), instead of just printing in the screen.
Thanks in advance.

[daniel klein]

It might be helpful to explain where blp is from. Without knowledge of what it is or does, you can return a scalar from Mata to Stata with

st_numscalar("name_for_stata_scalar", name_of_mata_scalar)

Note that name_for_stata_scalar will displace any existing Stata scalar with that name.

Best
Daniel

[Charlie Joyez]

It seems to me that the OP rather wants to return a variable into Stata than a scalar, but I might be wrong.
In that case the couple of lines below would do the trick:

mata st_addvar(type,newvarname)
mata st_store(.,newvarname,matamatrixname)

Where newvarname is the name of the variable you create (on Stata .dta), you use the same varname in the two lines. The matamatrixname is the original Mata (column) vector you want to store in Stata .dta
Note that the mata codes are not necessary if you already work under mata environment (but useful on a common Stata do-file).

See the example below

clear
set obs 3 /*your Stata .dta file should have at least the same number of observation than rows in your Mata vector*/

mata A=(1\2\3) /*creating a column vector in mata*/
mata A
mata st_addvar("double","var_A")
mata st_store(.,"var_A",A)
*Your variable is created in your .dta file!

Best,
Charlie

Edit : Actually, I figured this is quite odd if you work mainly on a Stata environment to create your Stata variable through a Mata command (first mata command), so please also consider another possibility (although the previous code is fine, this is why I keep the two versions):

clear
set obs 3
gen var_A=.

mata A=(1\2\3)
mata A
mata st_store(.,"var_A",A)

This way, you first create an empty variable in Stata, and then fill it the column vector stored in Mata.

[Bernardo Coelho]

Thanks Daniel and Charlie for your reply.
I tried what Charlie suggested but was not able to get the result I wanted.
I'll try to be more specific. Inside blp function, in the blp.ado, there is a function called eval_blp which is preceded by the command void. The variable I want is generated by this function.
I tried the st_addvar at the end of this eval_blp function, as I wrote in red. The eval_blp command is copied bellow.
I don't know if the fact that this vector is generated not by the original blp command, but by this command within a command makes difference.
Also, I didn't copy the whole blp command as it would be too big,
Thanks again!


void eval_blp(todo,theta,s,X,X_S,Z,W,noisily,Q,g,H)
{
//constructs GMM-objective function

real matrix X_S_t,E_t
real vector s_hat_t,ST,FI,w_t,delta
real scalar i,j,r,c,t

external D,beta,A,PI_HAT,R,tol_in,tol_out,D_A_delta,F,mkt_r ows,u
external N,V,T,mkt_rows_D,R,e_R,e_D,mv,w,K_1,K_S,K_D,K_PI,K _2

A=diag(exp(theta[1::K_S])) //trial values of random parameters in levels
_diag(A,editvalue(diagonal(A),0,1E-14))


//model includes demographic variables
if(!missing(D)){
e_D=J(K_D,1,1)
PI_VEC=theta[K_S+1::cols(theta)] //trial values of demographic parameters
PI_HAT=st_matrix("PI")

//create matrix of demographic parameters
c=0
for(i=1;i<=rows(PI_HAT);i++){
for(j=1;j<=cols(PI_HAT);j++){
if(PI_HAT[i,j]!=.){
c++
PI_HAT[i,j]=PI_VEC[c]
}
}
}
_editmissing(PI_HAT,0)
}

D_A_delta=J(N,K_2,.) //JT X K_2 matrix of derivatives of mean utility wrt heterogeneity parameters
ST=mkt_rows[,1] //market start rows
FI=mkt_rows[,2] //market finish rows

if(noisily){
displayas("txt")
printf("{txt:Contraction mapping by market}\n") //display iteration log
}

//computations by market
mv=0
for(t=1;t<=T;t++){
if(noisily){
t
}
tol=tol_in
X_S_t=panelsubmatrix(X_S,t,mkt_rows) //market observations for variables with random parameters
panelsubview(s_t,s,t,mkt_rows) //product shares
w_t=panelsubmatrix(w,t,mkt_rows)
E_t=randomEffects(X_S_t,t) //individual draws v_ijt where d_ijt=d_jt+v_ijt(theta2)
e_J=J(1,rows(w_t),1)
s_hat_t=((w_t:*E_t):/(e_J*(w_t:*E_t):+1))*e_R*(1/R) //predicted market shares
r=0

//contraction mapping
while (mreldif(ln(s_hat_t),ln(s_t))>tol){
w_t=w_t:*(s_t:/s_hat_t) //exponential of mean utiltiies
s_hat_t=((w_t:*E_t):/(e_J*(w_t:*E_t):+1))*e_R*(1/R)
++r
tol=tol_in*10^floor(r/1000) //reduces the tolerance by 10E-1 every 1000 iterations
}

w[|ST[t],.\FI[t],.|]=w_t //populate JT vector

D_A_delta[|ST[t],.\FI[t],.|]=partialDelta(w_t,E_t,X_S_t,t) //analytical derivatives
}

mata A = s_hat_t
mata st_addvar("double","s_hat_t")
mata st_store(.,"s_hat_t",A)

delta=log(w) //mean utilities
beta=invsym(X'Z*W*Z'X)*(X'Z*W*Z'delta) //estimator of parameters in mean utility equation
u=delta-X*beta //JT X 1 vector of residuals

//compute GMM objective function
if(mv){
Q=10E+10 //if not defined
}
else{
Q=editmissing((Z'u)'W*(Z'u),10E+10) //objective function
}

if(todo) {
g=(D_A_delta'*Z*W*Z'*u)' //analytical derivatives

}



if(noisily){
displayas("txt")
printf("{txt:Estimated standard deviations of random coefficients}\n")
A
if(!missing(D)){
printf("{txt:Estimated impacts of demographic variables (rows=coefficient equations, cols=variables)}\n")
PI_HAT
}
printf("{txt:Gradient of objective function}\n")
g
}

_diag(F=J(K_1+K_2,K_1+K_2,0),(J(1,K_1,1),exp(theta[1::K_S]),J(1,K_PI,1))) //delta method to compute var(theta)

}

[Sergio Correia]

As you are already in mata, you don't need to prefix the lines in red with "mata"

Also, it's useful to wrap the code in the

code

block, which makes it monospace and easier to read (it's the # button on the text bar above the editor)

[daniel klein]

I tried what Charlie suggested but was not able to get the result I wanted.

Please review the FAQ on reporting a problem. What do you mean by this statement? Do you get an error message? If so, which? If not, what exactly happens?

I would contact the author of blp and ask for help as there are just too many things that can go wrong here. Creating the variable in Mata will fail if blp preserves and later restores the original data. So you should probably create the variable before calling blp as Charlie suggested in his edited post. Further blp marks the observations to be used but this information is not readily available in the function you want to modify. Thus, filling the values with st_strore() without providing the information which observations are to be filled might give you wring results without you ever realizing it. There might be more complications ahead and I lack the time to get deeper into this.

Best
Daniel

[Charlie Joyez]

Also, you are not forced to use the A matrix as in my previous example.
This can cause troubles if you had already a mata matrix named A, and you still need it afterwards (because it would be erased), and it seems to be the case here, since I see those lines in the beginning of your code:

A=diag(exp(theta[1::K_S])) //trial values of random parameters in levels
_diag(A,editvalue(diagonal(A),0,1E-14))

So instead of the red lines you've wrote, try:

st_addvar("double","s_hat_t")
st_store(.,"s_hat_t",s_hat_t)

This way you directly use the vector s_hat_t defined in mata, without creating a new matrix A (and also removing mata prefix as suggested by Sergio)
However, are you sure this is a column vector? If not, transpose it in another matrix, but make sure the name isn't used already.

Concerning the deeper troubles in blp, I won't be useful and I suggest you to follow Daniel advice.
Please not that Daniel is also right when asking the precise error messages Stata displays.

Best,
Charlie

[Bernardo Coelho]

Thank you all for the replies, and sorry again for my inexperience in reporting the problem.
When I tried the code that Charlie had suggested, I got no error message. The command completed without any problem, but without saving the variable I wanted.
I tried now the:

st_addvar("double","s_hat_t")
st_store(.,"s_hat_t",s_hat_t)

But again, no error message and no variable created.
Thanks again,
Bernardo

[Charlie Joyez]

Dear Bernardo,

It is hard to believe that the command didn't work, and no error message displays, even in a quiet mode as I assume you are.
Make sure the command run (for example that the if condition under which the command runs, is verified), and the same for your initial vector s_hat_t in mata.

Could you please report the ouput of (after your code ran)

(mata) s_hat_t

It will tell wether the initial vector is created (only write the mata prefix if you're on Stata)

Best,
Charlie

[Bernardo Coelho]

Charlie,
I'll report below the output of the command. This command actually takes a long time to run with my full database (4 days), so I'm running a ligher version with only one market, witch is why convergence is not achieved, but right now i'm only interested in checking if I can retrieve s_hat_t.

blp s w3 w4, stochastic(x1, w3, cons) endog(x1=z1 z1w1 z1w2 z1w3 w12 w22 w1w2 w1w3 w2w3 z1w4 w1w4 w3w4 w2w4) markets(mun) draws(1000)
Iteration 0:   f(p) =  3.081e-33  
Iteration 1:   f(p) =          0  (not concave)
Iteration 2:   f(p) =          0  (not concave)
Iteration 49:  f(p) =          0  (not concave)
Iteration 50:  f(p) =          0  (not concave)
convergence not achieved

GMM estimator of BLP-model

GMM weight matrix: unadjusted               Number of obs            =  1
                                            Number of markets        =  1
                                            Number of Halton draws   =  1000
------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
Mean utility |
        cons |          0  (omitted)
          w3 |   .2764309          .        .       .            .           .
          w4 |          0  (omitted)
          x1 |          0  (omitted)
-------------+----------------------------------------------------------------
x1           |
          SD |         .5          .        .       .            .           .
-------------+----------------------------------------------------------------
w3           |
          SD |    .499922          .        .       .            .           .
-------------+----------------------------------------------------------------
cons         |
          SD |   .4999653          .        .       .            .           .
------------------------------------------------------------------------------


mata s_hat_t

 <istmt>:  3499  s_hat_t not found

Thanks again for all your help.
Best,
Bernardo

[Charlie Joyez]

Ok, as I though, your initial vector s_hat_t is not defined in mata, so you can't "convert" it in Stata variable (or more precisely, you can give its values to a Stata variable).
This is why no variable is created, but not error message displays.

We don't have your full code, but check the if conditions before introducing the s_hat_t, they might not all be verified.
Make sure your Mata code sets this vector.

Best,
Charlie

[Bernardo Coelho]

Dear Charlie,
Your explanation makes sense, but I thought the code was already defining s_hat_t.
The first time is appears is in:

 real vector s_hat_t,ST,FI,w_t,delta

Then it is defined by market by the iteration:

//computations by market
 mv=0
 for(t=1;t<=T;t++){        
         if(noisily){
        t
    }
    tol=tol_in
    X_S_t=panelsubmatrix(X_S,t,mkt_rows)            //market observations for variables with random parameters
    panelsubview(s_t,s,t,mkt_rows)                //product shares
    w_t=panelsubmatrix(w,t,mkt_rows)            
    E_t=randomEffects(X_S_t,t)                    //individual draws v_ijt where d_ijt=d_jt+v_ijt(theta2)
    e_J=J(1,rows(w_t),1)
    s_hat_t=((w_t:*E_t):/(e_J*(w_t:*E_t):+1))*e_R*(1/R)    //predicted market shares
    r=0

    //contraction mapping
    while (mreldif(ln(s_hat_t),ln(s_t))>tol){     
        w_t=w_t:*(s_t:/s_hat_t)                            //exponential of mean utiltiies
        s_hat_t=((w_t:*E_t):/(e_J*(w_t:*E_t):+1))*e_R*(1/R)
        ++r
        tol=tol_in*10^floor(r/1000)                        //reduces the tolerance by 10E-1 every 1000 iterations
    }
    
    w[|ST[t],.\FI[t],.|]=w_t                            //populate JT vector 
        
    D_A_delta[|ST[t],.\FI[t],.|]=partialDelta(w_t,E_t,X_S_t,t)        //analytical derivatives 
    
}

Aren't the parts in red the vector being defined in mata?
What should I code to make it be defined?

I cannot thank you enough, you are being a huge help.
Best,
Bernardo

[Xiaojin Wang]

Hi Bernardo,

I used -blp in Stata, and it took a long time to run. If you know Matlab, there're codes available for the blp model.

BLP code (Nevo original)
http://faculty.wcas.northwestern.edu...rc_dc_code.htm

BLP code 2 (revised by Hall, Berkeley)
http://eml.berkeley.edu/~bhhall/e220c/rc_dc_code.htm

Hope it can help.

Xiaojin