optimize: set skewness to zero

[Demetris Christodoulou]

Dear Statalisters,

I wish to write a simple optimisation problem to find the optimal fractional power for which the skewness of the targeted variable will be minimized (go to zero). I am aware of the commands lnskew0 and bcskew0. My question, although related to the idea of Tukey's ladder of powers, is more fundamental and related to the optimize code. Here is my code and output:

mata :
: void skew0(todo,p,y,g,H)
> {
> real vector X, Xp, d
> real scalar s
>
> X = st_data(.,"var1")
> Xp = X :^p
>
> d = Xp :-mean(Xp)
> y = mean(d:^3) / (mean(d:^2)):^(3/2)
> }
note: argument todo unused
note: argument g unused
note: argument H unused
note: variable s unused

: S = optimize_init()
: optimize_init_which(S,"min")
: optimize_init_evaluator(S, &skew0())
: optimize_init_technique(S,"nr")
: optimize_init_params(S,1)
: p = optimize(S)
Iteration 0: f(p) = 3.2341254
Iteration 1: f(p) = 3.0366119
Iteration 2: f(p) = .92755502 (not concave)
Iteration 3: f(p) = -.45619791 (not concave)
Iteration 4: f(p) = -.46533213
Iteration 5: f(p) = -.47138751 (not concave)
Iteration 6: f(p) = -.47154085 (not concave)
Iteration 7: f(p) = -.4715562
Iteration 8: f(p) = -.47155913 (not concave)
Iteration 9: f(p) = -.47155979

: printf("%12.8f",p) // iteration 0 gives the initial skewness of variable hin
0.00000004
: end

I know that the calculation of skewness works because iteration 0 gives the correct initial skewness of variable "var1" shown above at 3.2341254. However, I also know that the optimal value of p at 0.00000004 is incorrect as the skewness flips into negative and converges to -.47155979. I also know from other work that this power should be about 0.10.

Could someone more experienced that myself look at this code and help me understand what I am doing wrong.

many thanks in advance, Demetris Christodoulou

[German Rodriguez]

I would try modifying the criterion function to return the square of skewness. That would make skewness itself as close to zero as possible. You could modify the last line in skew0 to read

s = mean(d:^3) / (mean(d:^2)):^(3/2) 
y = s^2

[Demetris Christodoulou]

Excellent feedback German! This simple fix worked charms.