Stata Iteration Loop

[Irem Erten]

I have trouble completing a loop in which I set a parameter K (which can take "any" value), define a long series of functions, and have to make the output equal to the initial K until the iteration converges. how could I do this? I would like to give a simplified example as below.

gen K=0

forvalues i=0(0.0001)100 {
replace K=`i'
gen s=K^8/100-5+K^2

}

How can I adjust this code so that it eventually gives me which value of K makes s (the output of the entire code) equal to K? Thanks so much for any feedback! Really stuck here ....

[Clyde Schechter]

It looks like you are trying to solve the equation k = k^8/100 - 5 + k^2 = F(k) by a grid search.

At best this is a very inefficient way to solve this equation. The code uses variables that serve no purpose as you are really just working with single numbers, not series of observations.

Here's a fairly simple way. The value of k you seek also satisfies the equation f(k) = k^8/100 -5 + k^2 - k = 0. You can get a fairly fast convergence on a solution by starting with a very low value of k (your code suggested 0) and a very high one (your code suggests 100) and then you can iteratively pick the mid point of the interval between them. Evaluate f(mid point). If it is > 0, midpoint becomes the new upper value of k, if it is < 0, the midpoint becomes the new lower value of k. Keep going until the distance between the upper and lower values is within your desired precision (0.0001) and you will have it.

local klow = 0
local khigh = 100

foreach x in low high {
    local d`x' = `k`x''^8/100 - 5 + `k`x''^2 - `k`x''
}

display `dlow', `dhigh'
assert sign(`dhigh') == 1 & sign(`dlow') == - 1


while `khigh' - `klow' > 0.0001 {
    local k = (`khigh'+`klow')/2
    local d = `k'^8/100 -5 + `k'^2 - `k'
    if `d' > 0 {
        local khigh = `k'
    }
    else if `d' < 0 {
        local klow = `k'
    }
}
display `k', `k'^8/100 - 5 + `k'^2

The final display statement shows both k and F(k) so that you can see that they are nearly equal.

A more computationally efficient approach would be to use Newton's method, which usually converges much faster. But if you have a lot of functions to do this with, that would require you to find and program the derivative of each one, which probably would erase the efficiency gains. If memory serves, this method of halving the interval is due to Bolzano and Weierstrass.

Added:
A couple of things to think about here.

1. The function you show is an even-degree polynomial. In general, even-degree polynomials may not have any real roots. So beware of that.

2. The starting values of 0 and 100 seem to work here. But if the solution (assuming one even exists) does not lie between them, these iterations will just converge to one of those two numbers and leave you with no answer. So you need to have starting points that a) straddle an actual solution, and, b) have f(k) values with opposite signs. Problems can also arise if your function has two or more roots between the starting values. So the choice of starting values here is quite important.

3, If the function is such that at the initial high and low values the signs are the opposite of this example, then the code would need to be revised to switch the conditions under which khigh is reset to k and when k low is.

4. The loop terminating condition used here will cause the computations to end when klow and khigh are within 0.0001 of each other, which, assuming we are on the road to convergence, means that k will be within 0.0001 of the solution. However, if F(k) is steep near the actual root, you might find that F(k) is appreciably different from k. In the code above, our solution has k equaling F(k) to just 2 decimal places. Depending on why you are doing this in the first place, it may be more important to have F(k) very close to k then to have k within 0.0001 of the true value. So, you may have to impose a more stringent loop terminating condition to get that. For example, to get F(k) equal to k within 7 decimal places, you need to set the termination condition at 0.00000001.

[Nick Cox]

Cross-posted at http://stackoverflow.com/questions/4...iteration-loop

Please note our cross-posting policy, which is that you tell us about it. See the FAQ Advice.

[Scott Merryman]

You could use Mata's polyroot function which returns the real root of 2.0317481 and -1.66109948:

. mata: polyroots((-5, -1,1,0,0,0,0,0,.01))
                               1                           2                           3                           4                           5
    +---------------------------------------------------------------------------------------------------------------------------------------------
  1 |  -.076014341 - 2.39932716i   -.076014341 + 2.39932716i    1.70868001 - 1.49223156i    1.70868001 + 1.49223156i   -1.81798997 - 1.30021272i
    +---------------------------------------------------------------------------------------------------------------------------------------------
                               6                           7                           8
     -------------------------------------------------------------------------------------+
  1    -1.81798997 + 1.30021272i                   2.0317481                 -1.66109948  |
     -------------------------------------------------------------------------------------+