[Attaullah Shah]

Julia Schmidt I understand that you cannot share actual data, but can you share a dummy data set that has all the variables where the code and be seen working on them.

[Ingrid Qiu]

Hua Peng (StataCorp) I just noticed that I was replied. Thank you for your explanation. I made a mistake, thinking tolerance is the initial value. Now I know better about this function.

[Julia Schmidt]

Attaullah Shah unfortunately, I can't share the data but this is the code:

 generate icc_ct_=.
  mata
  z=J(1,1,.)
  st_view(z,., "icc_ct_ price bv_t0 feps_t1 feps_t2 feps_t3 feps_t4 feps_t5 growth payoutratio")
  function y(x,p,b0,f1,f2,f3,f4,f5,g,dv) {
      tf = f1+f2+f3+f4
      t1 = (f1-x*b0)/(1+x)
      t2 = (f2-x*(b0+f1*(1-dv)))/(1+x)^2
      t3 = (f3-x*(b0+(1-dv)*(f1+f2)))/(1+x)^3
      t4 = (f4-x*(b0+(1-dv)*(f1+f2+f3)))/(1+x)^4
      t5 = (f5-x*(b0+(1-dv)*tf))/(1+x)^5
      tv = ((f5-x*(b0+(1-dv)*tf))*(1+g))/((x-g)*(1+x)^5)
      return(-p+b0+t1+t2+t3+t4+t5+tv)
  }
  for (i=1;i<=rows(z);i++) {
      r=mm_root(icc_ct_=.,&y(),smallestdouble(),z[i,9]-epsilon(1),1e-9,1000,z[i,2],z[i,3],z[i,4],z[i,5],z[i,6],z[i,7],z[i,8],z[i,9],z[i,10])
      z[i,1]=icc_ct_
  }
  end
  
  
  generate icc_ct=.
  mata
  z=J(1,1,.)
  st_view(z,., "icc_ct price bv_t0 feps_t1 feps_t2 feps_t3 feps_t4 feps_t5 growth payoutratio")
  function s(x,p,b0,f1,f2,f3,f4,f5,g,dv) {
      tf = f1+f2+f3+f4
      t1 = (f1-x*b0)/(1+x)
      t2 = (f2-x*(b0+f1*(1-dv)))/(1+x)^2
      t3 = (f3-x*(b0+(1-dv)*(f1+f2)))/(1+x)^3
      t4 = (f4-x*(b0+(1-dv)*(f1+f2+f3)))/(1+x)^4
      t5 = (f5-x*(b0+(1-dv)*tf))/(1+x)^5
      tv = ((f5-x*(b0+(1-dv)*tf))*(1+g))/((x-g)*(1+x)^5)
      return(-p+b0+t1+t2+t3+t4+t5+tv)
  }
  for (i=1;i<=rows(z);i++) {
      r=mm_root(icc_ct=.,&s(),z[i,9]+epsilon(1),1-epsilon(1),1e-9,1000,z[i,2],z[i,3],z[i,4],z[i,5],z[i,6],z[i,7],z[i,8],z[i,9],z[i,10])
      z[i,1]=icc_ct
  }
  end

And then I merged icc_ct_ and icc_ct. I still have to verify, whether I did it correctly, so I better don't post that part of the code

[Attaullah Shah]

@ Julia Schmidt Can you please share an example data from your data and the code that worked for you.

[Julia Schmidt]

Hua Peng (StataCorp) Thank you so much!! Now it's working

[Hua Peng (StataCorp)]

The problem is the functions eq1() and eq2() are not continues in the domain you specified. Note Y:<y is discontinues at y. The following shows eq1() changes around 0.02 when the input is changed of 1e-15 at Y[38, 1]. Actually, function eq1() is discontinues at any observation Y where variable val is not 0.

. mata:
------------------------------------------------- mata (type end to exit) -----
: x = Y[38, 1]

: x
  6.420057774

: r1 = eq1(x, Y, val, P, t)

: r1
  -.2658067825

:
: r2 = eq1(x+(1e-15), Y, val, P, t)

: r2
  -.2411065604

: end
-------------------------------------------------------------------------------

mm_root() implements Brent's method https://en.wikipedia.org/wiki/Brent's_method which does assume the underlying function to be continues. That said, the algorithm is robust enough and can often deal with discontinues functions which behave reasonably. But the results need be carefully examined. In your above code, the first call to mm_root() is using a tolerance of 4.5 in a range 0-9 which seems way too large. Change it to 1e-5 gives the following:

: rc = mm_root(x1=.,&eq1(),0,9,1e-5,1000,Y,val,P,t)

: x1
  6.645090348

: rc
  0

:
: //Plot it back
: y = eq1(x1, Y, val, P, t)

: y
  -.0004896063

Use the same parameters but change the search interval to 6-9 give the same result:

: rc = mm_root(x1=.,&eq1(),6,9,1e-5,1000,Y,val,P,t)

: x1
  6.645088397

: rc
  0

:
: //Plot it back
: y = eq1(x1, Y, val, P, t)

: y
  -.0004896063

mm_root() believes it finds the root since rc is 0 and the result seems to reasonable when plotting back.

For eq2(),

: rc = mm_root(x1=.,&eq2(),0,9,1e-5,1000,Y,val,P,t)

: x1
  6.657811668

: rc
  0

:
: //Plot it back
: y = eq2(x1, Y, val, P, t)

: y
  -.0102072419

Lowering the tolerance and increasing the iteration do not improve the result:

: rc = mm_root(x1=.,&eq2(),0,9,1e-15,10000,Y,val,P,t)

: x1
  6.657814026

: rc
  0

:
: //Plot it back
: y = eq2(x1, Y, val, P, t)

: y
  -.0102072419

mm_root() believes it finds the root since rc is 0, but the results from plotting back raises concern which may be caused by the discontinuity of the underlying function.

[Ingrid Qiu]

Hua Peng (StataCorp) Julia Schmidt Thank you for your interest.

use data

mata:

st_view(Y=.,.,"Y")
st_view(val=.,.,"val")

P = 0.9676741647
t = 0.5
function eq1(y,Y,val,P,t) return(sum((Y:<y):*val)/sum(val) - (t-1+P)/P)
rc = mm_root(x1=.,&eq1(),0,9,4.5,1000,Y,val,P,t)

function eq2(y,Y,val,P,t) return(sum((Y:<y):*val)/sum(val) - t/P)
rc = mm_root(x2=.,&eq2(),0,9,0,1000,Y,val,P,t)

The above code gives the correct solution. I used

  quant = mm_quantile(Y,1,t)

as the initial value, but this will yield a wrong solution. Also, if I change the interval to 6-9, mm_root() will also return the wrong solution. I am not sure what makes this problem.

[Hua Peng (StataCorp)]

In your function, the last term is

tv = ((f5-x*(b0+(1-dv)*tf))*(1+g)/((x-g)*(1+x)^5)

which means the function is not defined when x = g, i.e., at the value of the growth variable. In Paul's version, the function is not defined at x=0. Hence you might need try the bounds from smallestdouble() to z[i,9]-epsilon(1); and z[i,9]+epsilon(1) to 1 if the value of growth variable is between 0 and 1.

[Julia Schmidt]

Hua Peng (StataCorp) The function is different because I try to estimate the ICC following Claus and Thomas, while Paul follows Gebhardt et al. The approaches differ in terms of the periods considered and how the terminal value is calculated. Since he was describing similar issues in his first post, I was hoping that someone could tell me, how to change my function to get the real results and not just the boundaries.
Ingrid Qiu what exactly is the initial value and how did you change it?
Thanks for your help

[Hua Peng (StataCorp)]

Julia Schmidt, your function is different from Paul's, it looks like the expansion of the same function but yours stops at 5th or 6th order, Paul's goes to 12th or 13th order, depending what the underlying function behaves given different coefficients "icc_ct price bv_t0 feps_t1 feps_t2 feps_t3 feps_t4 feps_t5 growth payoutratio", they may have very different roots. Again, this is just speculation without looking the coefficients.

Ingrid Qiu, I would be very interested if you could post the code of the problem you encountered.

[Ingrid Qiu]

Hi Stata Users,
I encountered a similar question when using mm_root(). I tried the method suggested above but it doesn't work for me. But after I changed the initial value, mm_root() returns the correct answers. I feel like mm_root() is very sensitive to the bound and the initial values. I am wondering whether there are something we need to pay attention to when we use this function. Thanks!

[Julia Schmidt]

Dear STATA Users,
I have a similar problem as Paul. I tried to estimate the ICC defined by Claus and Thomas (2001) as implemented by Li and Mohanram (2014). I used a similar approach as in Paul’s first posting and my results are mostly 0 or 1. There are some exceptions where STATA generates the correct value.
So I tried to change my STATA code according to your solution, Paul, but still, the results are not getting better.
Surprisingly, wenn I estimate ICC_GLS with your code, Paul, it generates the correct proxies.
I hope, someone can help me and finds the mistake in my code:

generate icc_ct=.
mata
z=J(1,1,.)
st_view(z,., "icc_ct price bv_t0 feps_t1 feps_t2 feps_t3 feps_t4 feps_t5 growth payoutratio")
function y(x,p,b0,f1,f2,f3,f4,f5,g,dv) {
    tf = f1+f2+f3+f4
    t1 = (f1-x*b0)/(1+x)
    t2 = (f2-x*(b0+f1*(1-dv)))/(1+x)^2
    t3 = (f3-x*(b0+(1-dv)*(f1+f2)))/(1+x)^3
    t4 = (f4-x*(b0+(1-dv)*(f1+f2+f3)))/(1+x)^4
    t5 = (f5-x*(b0+(1-dv)*tf))/(1+x)^5
    tv = ((f5-x*(b0+(1-dv)*tf))*(1+g)/((x-g)*(1+x)^5)
    return(-p+b0+t1+t2+t3+t4+t5+tv)
}
for (i=1;i<=rows(z);i++) {
    r=mm_root(icc_ct=.,&y(),smallestdouble(),1-epsilon(1),1e-9,1000,z[i,2],z[i,3],z[i,4],z[i,5],z[i,6],z[i,7],z[i,8],z[i,9],z[i,10])
    z[i,1]=icc_ct
}
end

Thanks a lot!
Julia


References:

• Claus, James and Jacob Thomas (2001): Equity Premia as Low as Three Percent? Evidence from Analysts’ Earnings Forecasts for Domestic and International Stock Markets, in: The Journal of Finance, Vol. 56, pp. 1629 – 1666.
• Li, Kevin K. and Partha Mohanram (2014): Evaluating Cross-Sectional Forecasting Models for Implied Cost of Capital, in: Review of Accounting Studies, Vol. 19, pp. 1152 – 1185.

[Paul Demere]

Thank you very much for the replies! On the advice of Hua Peng (StataCorp), I reached out to Ben Jann. He pointed out that the issue was with the bounds I was using, as the function is undefined at the bounds. By replacing the lower bound with smallestdouble() and the upper bound with 1-epsilon(1), I was able to run the code successfully. Included below is the code modified with the new bounds and made more readable at the request of Steve Samuels.

generate glsicc=.
mata
z=J(1,1,.)
st_view(z,., "glsicc p3prcc bve fore1 fore2 fore3 fore4 fore5 indmedroe dvrate")
function w(x,p,b0,f1,f2,f3,f4,f5,ii,dv) {
    tf = f1+f2+f3+f4
    fi = f5-ii
    t1 = (f1-x*b0)/(1+x)
    t2 = (f2-x*(b0+f1*dv))/(1+x)^2
    t3 = (f3-x*(b0+dv*(f1+f2)))/(1+x)^3
    t4 = (f4-x*(b0+dv*(f1+f2+f3)))/(1+x)^4
    t5 = (f5-x*(b0+dv*tf))/(1+x)^5
    t6 = ((f5-fi/7)-x*(b0+dv*(tf+f5)))/(1+x)^6
    t7 = ((f5-fi*2/7)-x*(b0+dv*(tf+2*f5-fi/7)))/(1+x)^7
    t8 = ((f5-fi*3/7)-x*(b0+dv*(tf+3*f5-fi*3/7)))/(1+x)^8
    t9 = ((f5-fi*4/7)-x*(b0+dv*(tf+4*f5-fi*6/7)))/(1+x)^9
    t10 = ((f5-fi*5/7)-x*(b0+dv*(tf+5*f5-fi*10/7)))/(1+x)^10
    t11 = ((f5-fi*6/7)-x*(b0+dv*(tf+6*f5-fi*15/7)))/(1+x)^11
    t12 = (ii-x*(b0+dv*(tf+7*f5-fi*3)))/((1+x)^12)
    tv = (ii-x*(b0+dv*(tf+7*f5-fi*3)))/(x*(1+x)^12)
    return(-p+b0+t1+t2+t3+t4+t5+t6+t7+t8+t9+t10+t11+t12+tv)
}
for (i=1;i<=rows(z);i++) {
    r=mm_root(glsicc=.,&w(),smallestdouble(),1-epsilon(1),1e-9,1000,z[i,2],z[i,3],z[i,4],z[i,5],z[i,6],z[i,7],z[i,8],z[i,9],z[i,10])
    z[i,1]=glsicc
}
end

[Hua Peng (StataCorp)]

Thanks Christophe Kolodziejczyk, I did not know -help moremata_source-.

[Christophe Kolodziejczyk]

You can look at the source code of the moremata package by typing help moremata_source. It is actually very instructive on Mata generally.

That being said you have very complicated function so it could be that the problem comes from here. An error in the coding of the function could mean that the problem has no solution. I would try first with a simpler function and with a dataset delivered with Stata. If you have problems there send the code and a log file of what has gone wrong.