Dear Statalisters,
For a dataset sorted on variable x, I wish to maximize the following function with respect to alpha. Here is the Latex equivalent:
\sum_i \sum_j [ \arctan(s_i / alpha) - \arctan(s_j / alpha) ]^2
where s_i = (\Delta y_i / \Delta x_i) and s_j = (\Delta y_j / \Delta x_j). The \Delta operator indicates first differences and s_i, s_j are consecutive slopes. Note how the optimization problem involves a double-summation operator over i=1,2,...,n and j=1,2,...,n-1, hence the need to create n*(n-1) pairs of s_i and s_j.
As an example, consider the sunspot dataset. Sort by x and calculate the first differences:
use http://www.stata-press.com/data/r14/sunspot.dta
rename spot y
rename time x
sort x
generate dy = y[_n]- y[_n-1]
generate dx = x[_n]- x[_n-1]
keep if dy<.
In Part 1 of the Mata code below, I take a view of dy and dx, and calculate all slopes dy/dx for every possible i,j pair where i!=j. I believe that this part works as expected.
In Part 2 of the Mata code, I decide to deal with double-summation by treating the i_th slopes as `panels' and the j_th slopes as `within panel' observations. I also reformulated the problem as Manhattan distances in a matrix form which should give the same result but it involves more complex code. To me, the `panel' approach seems to somewhat reduce the complexity of the problem, but please share your thoughts too.
My main concern is with the loop within the &eval function that is being optimized. I suspect that my literal translation of the optimization problem is inconsistent with Mata's optimize() engine. That is, for every i_th panel I calculate the colsum() of each panelsubmatrix() (i.e. the inner summation operator) and the loop accumulates for every i_th panel (i.e. the outer summation operator). Then, it is outside the loop that I specify that v should be optimized for the double-sum that is now reduced to just a scalar. This approach does not feel right to me, and also the end result is entirely unreasonable, but I am unable to see an alternative formulation.
Please share your wisdom. As this is part of a paper, all useful responses will be acknowledged accordingly.
many thanks, Demetris
mata:
mata clear
stata("sort " + "x")
st_view(dy=., ., "dy")
st_view(dx=., ., "dx")
// Part 1: get n*(n-1) pairs of slopes
id = J(rows(dy)*(rows(dy)-1),1,0)
si = J(rows(dy)*(rows(dy)-1),1,0)
sj = J(rows(dy)*(rows(dy)-1),1,0)
k = 0
for(i=1; i<=rows(dy); i++) {
for(j=1; j<=rows(dy); j++) {
if (i!=j) {
k=k+1
id[k] = i
si[k] = dy[i] / dx[i]
sj[k] = dy[j] / dx[j] ;
}
}
}
// Part 2: optimize
info = panelsetup(id, 1)
void eval(todo, a, info, si, sj, v, g, H)
{
sum_ij = 0
for (i=1; i<=rows(info); i++) {
sum_ij = sum_ij ///
+ colsum( panelsubmatrix( ( atan(si:/a) :- atan(sj:/a) ):^2, i, info ) )
}
v = sum_ij
}
S = optimize_init()
optimize_init_evaluator(S, &eval())
optimize_init_evaluatortype(S, "gf0")
optimize_init_which(S,"max")
optimize_init_params(S, 1)
optimize_init_argument(S, 1, info)
optimize_init_argument(S, 2, si)
optimize_init_argument(S, 3, sj)
a = optimize(S)
end
For a dataset sorted on variable x, I wish to maximize the following function with respect to alpha. Here is the Latex equivalent:
\sum_i \sum_j [ \arctan(s_i / alpha) - \arctan(s_j / alpha) ]^2
where s_i = (\Delta y_i / \Delta x_i) and s_j = (\Delta y_j / \Delta x_j). The \Delta operator indicates first differences and s_i, s_j are consecutive slopes. Note how the optimization problem involves a double-summation operator over i=1,2,...,n and j=1,2,...,n-1, hence the need to create n*(n-1) pairs of s_i and s_j.
As an example, consider the sunspot dataset. Sort by x and calculate the first differences:
use http://www.stata-press.com/data/r14/sunspot.dta
rename spot y
rename time x
sort x
generate dy = y[_n]- y[_n-1]
generate dx = x[_n]- x[_n-1]
keep if dy<.
In Part 1 of the Mata code below, I take a view of dy and dx, and calculate all slopes dy/dx for every possible i,j pair where i!=j. I believe that this part works as expected.
In Part 2 of the Mata code, I decide to deal with double-summation by treating the i_th slopes as `panels' and the j_th slopes as `within panel' observations. I also reformulated the problem as Manhattan distances in a matrix form which should give the same result but it involves more complex code. To me, the `panel' approach seems to somewhat reduce the complexity of the problem, but please share your thoughts too.
My main concern is with the loop within the &eval function that is being optimized. I suspect that my literal translation of the optimization problem is inconsistent with Mata's optimize() engine. That is, for every i_th panel I calculate the colsum() of each panelsubmatrix() (i.e. the inner summation operator) and the loop accumulates for every i_th panel (i.e. the outer summation operator). Then, it is outside the loop that I specify that v should be optimized for the double-sum that is now reduced to just a scalar. This approach does not feel right to me, and also the end result is entirely unreasonable, but I am unable to see an alternative formulation.
Please share your wisdom. As this is part of a paper, all useful responses will be acknowledged accordingly.
many thanks, Demetris
mata:
mata clear
stata("sort " + "x")
st_view(dy=., ., "dy")
st_view(dx=., ., "dx")
// Part 1: get n*(n-1) pairs of slopes
id = J(rows(dy)*(rows(dy)-1),1,0)
si = J(rows(dy)*(rows(dy)-1),1,0)
sj = J(rows(dy)*(rows(dy)-1),1,0)
k = 0
for(i=1; i<=rows(dy); i++) {
for(j=1; j<=rows(dy); j++) {
if (i!=j) {
k=k+1
id[k] = i
si[k] = dy[i] / dx[i]
sj[k] = dy[j] / dx[j] ;
}
}
}
// Part 2: optimize
info = panelsetup(id, 1)
void eval(todo, a, info, si, sj, v, g, H)
{
sum_ij = 0
for (i=1; i<=rows(info); i++) {
sum_ij = sum_ij ///
+ colsum( panelsubmatrix( ( atan(si:/a) :- atan(sj:/a) ):^2, i, info ) )
}
v = sum_ij
}
S = optimize_init()
optimize_init_evaluator(S, &eval())
optimize_init_evaluatortype(S, "gf0")
optimize_init_which(S,"max")
optimize_init_params(S, 1)
optimize_init_argument(S, 1, info)
optimize_init_argument(S, 2, si)
optimize_init_argument(S, 3, sj)
a = optimize(S)
end
