Calling Mata

[Runliang Li]

Hallo Everyone:

I have compiled a function in Mata and I can run it in the interactive mode (lasts about 30 minutes). This Mata function needs a series of matrices which I have already put into Mata using the "putmata" syntax. My problem is that I cannot execute the function if I call it in the "nlgrid" program. I tried to put the program and the mata into different files but that still did not work. The major part of my do-file is as follows. Thanks in advance for your comments and suggestions.

putmata estdate myid_bodqrt0 = (myid_bodqrt1 myid_bodqrt2 myid_bodqrt3 myid_bodqrt4) ageqrt0 = (ageqrt1 ageqrt2 ageqrt3 ageqrt4) Eqrt0 = (Eqrt1 Eqrt2 Eqrt3 Eqrt4) EBqrt0 = (EBqrt1 EBqrt2 EBqrt3 EBqrt4) denomqrt0 = (denomqrt1 denomqrt2 denomqrt3 denomqrt4) deposum omega myid cntynumb, replace

cap mata: mata drop eqexpbr()
version 14.0
mata:
void eqexpbr(real scalar k)
{

real matrix Bu, myid_bodqrt0, myid_blc, UU, estdate, ageqrt0, denomqrt0, eqall, C1, cntynumb, C2, S, eqlogmat, Eqrt0, omega, EB

Bu = uniqrows(myid_bodqrt0)
myid_blc = vec(Bu)
UU = J(rows(Bu),4,.)

for (q=1; q<=4; q++) {
for (i=1; i<=rows(estdate); i++) {
if (ageqrt0[i,q]>90) {
for (j=90; j<=ageqrt0[i,q]-1; j++) {
denomqrt0[i,q] = denomqrt0[i,q] + (ageqrt0[i,q]-j)^k
}
}

C1 = eqall[,1..4]
C2 = J(rows(C1),cols(C1),cntynumb[i])
S = C1-C2

Sprod = S[,1]:*S[,2]:*S[,3]:*S[,4]
pos = Sprod:==0 // pos = 1 if the countyfp is matched.
if (sum(pos) != 0) {
eqlogmat = select(eqall[,5..6],pos)
}
for(l=1; l<=rows(eqlogmat); l++) {
if (eqlogmat[l,2]>=estdate[i] & denomqrt0[i,q]>0 & Eqrt0[i,q]>=0) {
Eqrt0[i,q] = Eqrt0[i,q]+(eqlogmat[l,2]-estdate[i])^k/denomqrt0[i,q]*eqlogmat[l,1]
}
}
}

for (u=1; u<=rows(Bu); u++) {
UU[u,]=colsum(omega:*Eqrt0:*rowsum(myid_bodqrt0:==Bu[u,]))
}
}

EB = vec(UU)
stata("getmata EB, id(myid_blc) replace")
}
end


capture program drop nleqexp
program nleqexp
version 14.0
syntax [varlist(default=none)] if, at(name)
tokenize `varlist'
marksample touse
local n : word count `varlist'
local n = `n'-1 // no intercept in these cases
tempname a b k
local atn = colsof(`at') // includes expl var | learning-from-experience prelim coeff | intercept
local pn = `atn'-`n' // # of learning-from-experience parameters: b, c, k
scalar `k' = `at'[1,`n']

scalar `a' = 0
qui replace `1' = `a' `if' // token `1' is the depvar, here start replacing it with fitted value
forvalues i=2/`n' {
tempname ipar
scalar `ipar' = `at'[1,`i'-1] // retrieve parameter values
qui replace `1' = `1' + `ipar'*``i'' `if' // calculates function value
}
end


capture program drop nlgrid
program nlgrid, eclass
version 14.0
syntax varlist [if] [, KLOW(real 1) KSTEP(real 1) KHI(real 1)]
marksample touse
gettoken depvar explvar: varlist
tempname rss b minrss binit binitmod kmin v vinit
local parnl "k"

forvalues k = `klow'(`kstep')`khi' {
mata: eqexpbr(`k')
qui reg `varlist' if `touse', noconstant
scalar `rss' = e(rss)
matrix `b' = e(b)
matrix `v' = vecdiag(e(V))

if `k' == `klow' {
scalar `minrss' = `rss'
matrix `binit' = `b'
matrix `vinit' = `v'
scalar `kmin' = `k'
}
else if `rss' < `minrss' {
scalar `minrss' = `rss'
matrix `binit' = `b'
matrix `vinit' = `v'
scalar `kmin' = `k'
}

di as text "k " as result `k'
di as text "rss " as result `rss'
}
di as text "k " as result `kmin'
di as text "rss " as result `minrss'

local nvar = colsof(`binit')
local pvar : colfullnames(`binit') // extract variable names
tokenize `pvar' // this way variables that reg above dropped are also dropped
local explvar ""
qui replace `touse' = e(sample)

forvalues h = 1/`nvar' {
matrix `binitmod' = nullmat(`binitmod'),`binit'[1,`h']
local explvar "`explvar' ``h''"
}
matrix `binit' = `binitmod',`kmin'
nl eqexp @ `depvar' `explvar' if `touse', variables(`explvar') parameters(`explvar' `parnl') initial(`binit') // delta(4e-4)
end

eststo r1: xi: nlgrid EquityRatio rcon2170 ROA ImputedIntRateTT RERatio NonHazPortion Bankage i.YearReporting if (EquityRatio ~=.), klow(0.8) kstep(0.4) khi(4)

In addition, the related error messages are:

eststo r1: xi: nlgrid EquityRatio rcon2170 ROA ImputedIntRateTT RERatio NonHazPortion Bankage
> i.YearReporting if (EquityRatio ~=.), klow(0.8) kstep(0.4) khi(4)
i.YearReporting _IYearRepor_1994-2013(naturally coded; _IYearRepor_1994 omitted)

vector or matrix myid_blc not found
stata(): 3598 Stata returned error
eqexpbr(): - function returned error
<istmt>: - function returned error

[Sergio Correia]

It's hidden in the help file of putmata (at the end), but basically don't use putmata/getmata within an ado file or a program. Instead, use st_data() and st_store()

[Runliang Li]

Thank you Sergio. I can call the Mata function if I use st_data() and st_store(). By using the stata() syntax, I can also put data of different sources into Mata.

[Runliang Li]

Hallo Everyone:

I would like to continue this topic because I encounter a new problem when proceeding with my estimation and this problem is still related to calling Mata in the program. First at all, I state my ultimate purpose. Basically, I want to use the "nl" command to estimate a non-linear regression. The last variable, EB, is a relatively complicated function of a parameter k and data of other sources. I generate this variable by composing a Mata function "eqexpbr". In the initial value finding program "nlgrid", I try different values of k and decide to use the k that gives the smallest rss of the linear regression. In the non-linear program "nleqexp", I try to let Stata know that EB is a function of the parameter k so that Stata can carry out the iteration. However, I did not succeed. The following paragraphs are my programs. Thanks in advance for your advice and comments.

use $pathdata/HDV_Balanceandincome.dta, clear
drop if rssd9050 == 0
drop if missing(rssd9999)
sort rssd9050 DateReporting
qui ta rssd9050 // di `r(r)' shows that there are 628 banks.
tostring rssd9050, gen(rssd9050str)
gen myid_blc = rssd9050str+rssd9999 // concatenate i-id and t-id to just a single id so as to use getmata, id().
destring myid_blc, replace
label variable myid_blc "COMBINED ID, BLC STANDS FOR BALANCE SHEET."
qui tab YearReporting, gen(YearReporting)
gen EB = .
save $pathdata/HDV_Balanceandincome_mata.dta, replace

cap mata: mata drop eqexpbr()
version 14.0
mata:
void eqexpbr(real scalar k)
{

stata("use HDV_EQcountymatch_mata.dta, clear")
eqall = st_data(.,"countyfp_temp1 countyfp_temp2 countyfp_temp3 countyfp_temp4 magnum time")

stata("use HDV_BOD19942016_pureappend_mata.dta, clear")
estdate = st_data(.,"estdate")
myid_bodqrt0 = st_data(.,"myid_bodqrt1 myid_bodqrt2 myid_bodqrt3 myid_bodqrt4")
ageqrt0 = st_data(.,"ageqrt1 ageqrt2 ageqrt3 ageqrt4")
Eqrt0 = st_data(.,"Eqrt1 Eqrt2 Eqrt3 Eqrt4")
EBqrt0 = st_data(.,"EBqrt1 EBqrt2 EBqrt3 EBqrt4")
denomqrt0 = st_data(.,"denomqrt1 denomqrt2 denomqrt3 denomqrt4")
deposum = st_data(.,"deposum")
omega = st_data(.,"omega")
myid = st_data(.,"myid")
cntynumb = st_data(.,"cntynumb")

Bu = uniqrows(myid_bodqrt0)
myid_blc = vec(Bu)
UU = J(rows(Bu),4,.)

for (q=1; q<=4; q++) {
for (i=1; i<=rows(estdate); i++) {
if (ageqrt0[i,q]>90) {
for (j=90; j<=ageqrt0[i,q]-1; j++) {
denomqrt0[i,q] = denomqrt0[i,q] + (ageqrt0[i,q]-j)^k
}
}

C1 = eqall[,1..4]
C2 = J(rows(C1),cols(C1),cntynumb[i])
S = C1-C2

Sprod = S[,1]:*S[,2]:*S[,3]:*S[,4]
pos = Sprod:==0 // pos = 1 if the countyfp is matched.
if (sum(pos) != 0) {
eqlogmat = select(eqall[,5..6],pos)
}
for (l=1; l<=rows(eqlogmat); l++) {
if (eqlogmat[l,2]>=estdate[i] & denomqrt0[i,q]>0 & Eqrt0[i,q]>=0) {
Eqrt0[i,q] = Eqrt0[i,q]+(eqlogmat[l,2]-estdate[i])^k/denomqrt0[i,q]*eqlogmat[l,1]
}
}
}

for (u=1; u<=rows(Bu); u++) {
UU[u,]=colsum(omega:*Eqrt0:*rowsum(myid_bodqrt0:==Bu[u,]))
}
}

EB = vec(UU)

stata("use HDV_Balanceandincome_mata.dta, clear")
id_blc = st_data(.,"myid_blc")
eb_bod = J(rows(id_blc),1,.)
for (s=1; s<=rows(id_blc); s++) {
eb_bod[s] = colsum(EB:*(myid_blc:==id_blc[s]))
}
st_store(.,"myid_blc",id_blc)
st_store(.,"EB",eb_bod)
stata("save HDV_Balanceandincome_mata.dta, replace")
}
end


capture program drop nleqexp
program nleqexp, sortpreserve
version 14.0
syntax [varlist(default=none)] if, at(name)
tokenize `varlist'
//marksample touse
local n : word count `varlist' // this n is the number of all parameters because #para(X) = n-1, #nl = 1.
local n = `n'-1
tempname a b k

scalar `a' = 0
qui replace `1' = `a' `if' // token `1' is the depvar, here start replacing it with fitted value

local m = `n'-1

forvalues i=2/`m' {
tempname ipar
scalar `ipar' = `at'[1,`i'-1] // retrieve parameter values
qui replace `1' = `1' + `ipar'*``i'' `if' // calculates function value
}

scalar `k' = `at'[1,`n']
local kn = `k' // I can not implement mata: eqexpbr(`k') but have to create a local. Do you know why?
tempvar EBNL
preserve
mata: eqexpbr(`kn')
generate double `EBNL' = EB `if'
restore
qui replace `1' = `1' + `at'[1,`n'-1]*`EBNL'
end

capture program drop nlgrid
program nlgrid, eclass
version 14.0
syntax varlist [if] [, KLOW(real 1) KSTEP(real 1) KHI(real 1)]
tempname rss b minrss binit binitmod kmin v vinit
local parnl "k"

forvalues k = `klow'(`kstep')`khi' {
mata: eqexpbr(`k')
local varlistnew "`varlist' EB"
gettoken depvar explvar: varlistnew
marksample touse
qui reg `varlistnew' if `touse' , noconstant
scalar `rss' = e(rss)
matrix `b' = e(b)
matrix `v' = vecdiag(e(V))

if `k' == `klow' {
scalar `minrss' = `rss'
matrix `binit' = `b'
matrix `vinit' = `v'
scalar `kmin' = `k'
}
else if `rss' < `minrss' {
scalar `minrss' = `rss'
matrix `binit' = `b'
matrix `vinit' = `v'
scalar `kmin' = `k'
}

di as text "k " as result `k'
di as text "rss " as result `rss'
}
di as text "k " as result `kmin'
di as text "rss " as result `minrss'

local nvar = colsof(`binit')
local explvar ""
qui replace `touse' = e(sample)

forvalues h = 1/`nvar' {
matrix `binitmod' = nullmat(`binitmod'),`binit'[1,`h']
local explvar "`explvar' ``h''"
}
matrix `binit' = `binitmod',`kmin'
nl eqexp @ `depvar' `explvar' if `touse', variables(`explvar') parameters(`explvar' `parnl') initial(`binit') noconstant
end

// Main Body
eststo r1: nlgrid EquityRatio rcon2170 ROA ImputedIntRateTT RERatio NonHazPortion Bankage YearReporting1-YearReporting20 if (EquityRatio ~=.), klow(0.8) kstep(0.4) khi(2)

To be more detailed, the error message before using the "preserve" or "restore" is:

nleqexp returned 111
verify that nleqexp is a function evaluator program
could not restore sort order because variables were dropped

The error message after using the "preserve" or "restore" is:

nleqexp returned 111
verify that nleqexp is a function evaluator program

[Runliang Li]

Could anyone have a look at my problem?

[Sergio Correia]

Quick advice:

• Show us a minimum working example (mwe) of your issue, this is just too long to read
• Wrap the code in the CODE button on top of the textbox, it's way easier to read
• You CANNOT use stata("use..") within a mata function like that, because the nl evaluator likely saves stuff on the dataset which you then delete. And the error message is literally this: "could not restore sort order because variables were dropped"
• A soln might be to keep whatever thata you need to load as a mata matrix, which will also make the program much faster

S

[Runliang Li]

Hallo Sergio:

Can you tell me how to load data as a mata matrix? I have tried the following code:

mata:
eqall = st_data(.,"countyfp_temp1 countyfp_temp2 countyfp_temp3 countyfp_temp4 magnum time")
end
mata:
estdate = st_data(.,"estdate")
// other data omitted.
end

cap mata: mata drop eqexpbr()
version 14.0
mata:
function eqexpbr(scalar k)  
{
for (q=1; q<=4; q++)  {
   for (i=1; i<=rows(estdate); i++)  {
      if (ageqrt0[i,q]>90)  {
         for (j=90; j<=ageqrt0[i,q]-1; j++)  {
            denomqrt0[i,q] = denomqrt0[i,q] + (ageqrt0[i,q]-j)^k
         }
      }
 
   C1 = eqall[,1..4]
   C2 = J(rows(C1),cols(C1),cntynumb[i])
   S = C1-C2
   
   Sprod = S[,1]:*S[,2]:*S[,3]:*S[,4]
   pos = Sprod:==0 // pos = 1 if the countyfp is matched.
      if (sum(pos) != 0)  {
         eqlogmat = select(eqall[,5..6],pos)
      }
      for (l=1; l<=rows(eqlogmat); l++)  {
         if (eqlogmat[l,2]>=estdate[i] & denomqrt0[i,q]>0 & Eqrt0[i,q]>=0)  {  
            Eqrt0[i,q] = Eqrt0[i,q]+(eqlogmat[l,2]-estdate[i])^k/denomqrt0[i,q]*eqlogmat[l,1]
         } 
      }
   }
   
  for (u=1; u<=rows(Bu); u++)  {  
     UU[u,]=colsum(omega:*Eqrt0:*rowsum(myid_bodqrt0:==Bu[u,]))
   } 
}

EBm = vec(UU)

st_view(eb_bod=.,.,"EB")
for (s=1; s<=rows(id_blc); s++)  {
   eb_bod[s] = colsum(EBm:*(myid_blc:==id_blc[s]))
}
}
end

But when I try "mata: eqexpbr(2)", it does not go through and nothing is changed.

Regards,
RL

[Sergio Correia]

TBH not sure what is going on (although I appreciate the cleaner code syntax!)

You probably want to try it with a very simple case, including some debuging code that prints results or text on the window, so you know if the code was run.