-
Notifications
You must be signed in to change notification settings - Fork 8
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #5 from pcorralrodas/weighted_lnskew
Added lnskew_w option : Commands have been tested with new option
- Loading branch information
Showing
4 changed files
with
244 additions
and
8 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,214 @@ | ||
| *! version 1.0 23/02/2023 | ||
| * Original code from Stata - adjusted to allow for survey weights | ||
| cap prog drop lnskew0w | ||
| cap prog drop Lnskew0w | ||
| program define lnskew0w, rclass sort | ||
| version 6, missing | ||
| global S_1 . /* Gamma */ | ||
| global S_2 . /* Lower confidence limit */ | ||
| global S_3 . /* Upper confidence limit */ | ||
| global S_4 . /* skewness */ | ||
|
|
||
| syntax anything(name=eq equalok) [if] [in] [, Level(passthru) * weight(varname) ] | ||
| if "`level'" != "" { | ||
| local old `"`eq' `if' `in' , `options'"' | ||
| local 0 ", `level'" | ||
| syntax [, level(cilevel)] | ||
| local 0 `"`old'"' | ||
| local ci = `level'/100 | ||
| } | ||
| else { | ||
| local ci 0 | ||
| local level 95 | ||
| } | ||
| syntax newvarname =/exp [if] [in] [, /* | ||
| */ Delta(real 0) Zero(real 0) weight(varname)] | ||
|
|
||
| if ("`weight'"==""){ | ||
| tempvar weight | ||
| gen `weight' = 1 | ||
| } | ||
|
|
||
| local exp `=trim(`"`exp'"')' // rm white space, messes up output | ||
| tempvar x lnx | ||
| quietly { | ||
| gen double `lnx'=. | ||
| gen double `x' = `exp' `if' `in' | ||
| summ `x' [aw=`weight'], detail | ||
| local obs = r(N) | ||
| if `obs' < 3 { noisily error 2001 } | ||
| local skew = r(skewness) | ||
| local min = r(min) | ||
| local max = r(max) | ||
| if `min' == `max' { | ||
| noisily di in red "no variance" | ||
| exit 409 | ||
| } | ||
| local min0 = r(min) | ||
| local range = r(max)-r(min) | ||
| /* | ||
| Standardize data so that min=0, max=1. | ||
| */ | ||
| replace `x' = (`x'-`min0')/`range' | ||
| local min = 0 | ||
| local max = 1 | ||
| /* | ||
| Start of Gstar. | ||
| */ | ||
| local xmed = (r(p50)-`min0')/`range' | ||
| local term = `min'+`max'-2.0*`xmed' | ||
| local gamma = -4 | ||
| if `term'>0 { | ||
| local g = `min'-(`min'-`xmed')^2/`term' | ||
| } | ||
| else if `term'<0 { | ||
| local g =-`max'+(`max'-`xmed')^2/`term' | ||
| } | ||
| else { | ||
| local g -4 | ||
| } | ||
| local m1 = 1 | ||
| if `skew'<0 { | ||
| local min -1 | ||
| local max 0 | ||
| replace `x' = -`x' | ||
| local minus "-" | ||
| local m1 -1 | ||
| } | ||
| if (`g'>-4) & (`term'<0)==(`skew'<0) { local gamma = `g' } | ||
| /* | ||
| End of Gstar. | ||
| */ | ||
| local delta = cond(`delta'<=0, .02, `delta') | ||
| local zero = cond(`zero' <=0, .001, `zero') | ||
| local target 0 | ||
| Lnskew0w `x' `gamma' `min' `delta' `zero' `target' `lnx' `weight' | ||
| local gamma = r(gamma) /* gamma in transformed units */ | ||
| local skewg = r(skewness) | ||
| /* | ||
| Find confidence interval for gamma in transformed units | ||
| if n>=8. | ||
| */ | ||
| if (`obs'>7) & (`ci'>0) { | ||
| local z = invnorm(0.5+0.5*`ci') | ||
| _crczsku `z' `obs' 1 | ||
| local target = r(target) | ||
| /* | ||
| Lower confidence limit (if it exists). | ||
| If not, it is taken as missing (minus infinity). | ||
| */ | ||
| if abs(`skew')<`target' { local gammal . } | ||
| else { | ||
| Lnskew0w `x' `gamma' `min' /* | ||
| */ `delta' `zero' `target' `lnx' `weight' | ||
| local gammal = `m1'*`min0'+r(gamma)*`range' | ||
| } | ||
| local target = -`target' | ||
| Lnskew0w `x' `gamma' `min' `delta' /* | ||
| */ `zero' `target' `lnx' `weight' | ||
| local gammah = `m1'*`min0'+r(gamma)*`range' | ||
| } | ||
| else local ci 0 | ||
| replace `lnx'=`lnx'+log(`range') | ||
| replace `lnx'=log(`x'-`gamma')+log(`range') | ||
| local gamma = `m1'*`min0'+`gamma'*`range' | ||
| } | ||
| local cil `=string(`level')' | ||
| local cil `=length("`cil'")' | ||
| local spaces "" | ||
| if `cil' == 2 { | ||
| local spaces " " | ||
| } | ||
| else if `cil' == 4 { | ||
| local spaces " " | ||
| } | ||
| di _n in smcl in gr _col(8) "Transform" _col(18) "{c |}" /* | ||
| */ _col(29) "k" _col(31) "`spaces'" /* | ||
| */ `"[`=strsubdp("`level'")'% conf. interval] Skewness"' _n /* | ||
| */ "{hline 17}{c +}{hline 50}" | ||
| local lhs = "ln(`minus'"+bsubstr("`exp'",1,8)+"-k)" | ||
| local l=16-length("`lhs'") | ||
| di in smcl in gr _skip(`l') "`lhs' {c |} " /* | ||
| */ in ye %9.0g `gamma' " " _c | ||
| if `ci' { | ||
| ret scalar lb = `gammal' | ||
| ret scalar ub = `gammah' | ||
| global S_2 `gammal' | ||
| global S_3 `gammah' | ||
| if `gammal'>=. { | ||
| di "-infinity" _c | ||
| } | ||
| else di %9.0g `gammal' _c | ||
| di " " %9.0g `gammah' _c | ||
| } | ||
| else { | ||
| di in gr " (not calculated) " _c | ||
| } | ||
| ret scalar skewness = `skewg' | ||
| global S_4 `skewg' | ||
| di " " %9.0g return(skewness) | ||
| ret scalar gamma = `gamma' | ||
| global S_1 `gamma' | ||
| local gamma : display string(`gamma', "%9.0g") | ||
| local gamma = trim(`"`gamma'"') | ||
| if `gamma' > 0 { | ||
| local label `"ln(`minus'`exp'-`gamma')"' | ||
| } | ||
| else { | ||
| local gamma = -`gamma' | ||
| local label `"ln(`minus'`exp'+`gamma')"' | ||
| } | ||
| gen `typlist' `varlist' = `lnx' | ||
| label var `varlist' `"`label'"' | ||
| end | ||
|
|
||
|
|
||
| program define Lnskew0w, rclass | ||
| /* | ||
| Args: 1=_x, 2=gamma, 3=`min', 4=`delta', 5=`zero', 6=target skewness. | ||
| 7=_lnx variable | ||
| Returns gamma in r(gamma), skewness in r(skewness) | ||
| */ | ||
| args X | ||
| local gamma = `2' | ||
| local min = `3' | ||
| local delta = `4' | ||
| local zero = `5' | ||
| local target = `6' | ||
| local lnx "`7'" | ||
| local peso "`8'" | ||
|
|
||
| local iter 0 | ||
| local eta = -log(`min'-`gamma') | ||
| replace `lnx' = log(`X'-`gamma') | ||
| sort `lnx' /* for speed */ | ||
| summ `lnx' [aw=`peso'], detail | ||
| local f0 = r(skewness)-`target' | ||
| ret scalar skewness = r(skewness) | ||
| * di in red "Iter Gamma Skewness" | ||
| while (abs(`f0') > `zero') { | ||
| * noisily di in red `iter' " " `gamma' " " `f0' | ||
| local iter = `iter'+1 | ||
| replace `lnx' = log(`X'-`min'+exp(-`eta'-`delta')) | ||
| sum `lnx' [aw=`peso'], detail | ||
| local m = (r(skewness)-`target'-`f0')/`delta' | ||
| if `m' == 0 { | ||
| * noisily di in red /* | ||
| */ "(convergence problems, doubling the value of delta)" | ||
| local delta = `delta'*2 | ||
| } | ||
| else { | ||
| local eta = `eta'-`f0'/`m' | ||
| local gamma = `min'-exp(-`eta') | ||
| if (`gamma' > `min') { | ||
| local gamma = `min'-.0000001 | ||
| local eta = -log(`min'-`gamma') | ||
| } | ||
| replace `lnx' = log(`X'-`gamma') | ||
| sum `lnx' [aw=`peso'], detail | ||
| local f0 = r(skewness)-`target' | ||
| ret scalar skewness = r(skewness) | ||
| } | ||
| } | ||
| ret scalar gamma = `gamma' | ||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters