Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Extracting results from MATA

    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.

  • #2
    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

    Code:
    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

    Comment


    • #3
      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:
      Code:
      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

      Code:
      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):
      Code:
      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.
      Last edited by Charlie Joyez; 01 Nov 2016, 14:14.

      Comment


      • #4
        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)

        }

        Comment


        • #5
          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:
          code
          block, which makes it monospace and easier to read (it's the # button on the text bar above the editor)

          Comment


          • #6
            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
            Last edited by daniel klein; 02 Nov 2016, 00:11.

            Comment


            • #7
              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:
              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:
              Code:
              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

              Comment


              • #8
                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:
                Code:
                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

                Comment


                • #9
                  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)
                  Code:
                  (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

                  Comment


                  • #10
                    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.
                    Code:
                    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

                    Comment


                    • #11
                      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

                      Comment


                      • #12
                        Dear Charlie,
                        Your explanation makes sense, but I thought the code was already defining s_hat_t.
                        The first time is appears is in:
                        Code:
                         real vector s_hat_t,ST,FI,w_t,delta
                        Then it is defined by market by the iteration:

                        Code:
                        //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

                        Comment


                        • #13
                          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

                          Comment

                          Working...
                          X