Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: Re: Stuck at akdensity

From   Ramani Gunatilaka <>
Subject   st: Re: Stuck at akdensity
Date   Sun, 07 Dec 2003 01:41:49 +0000

Hi all,
Stuck again this weekend, I'm afraid. I am trying to plot akdensity distributions using Philippe Van Kerm's SSC download. 
My programme gets stuck within the akdensity routine, and I can't figure out why. 

This is the log output of set trace on immediately preceding error message.

   ------------------------------------------------------- end ipolate ---
      - if "`zero'"~="" {replace ynew = 0 if ynew==. & x~=. }
      = if "zero"~="" {replace ynew = 0 if ynew==. & x~=. }
      - }
      - qui keep if _stack==2
      - rename ynew `tmp'
      = rename ynew __00000K
      - drop y x
      - rename pos `tmppos'
      = rename pos __00000L
      - drop _stack
      - tempfile tmpsav
      - sort `tmppos'
      = sort __00000L
      - qui save `tmpsav' , replace
      = qui save C:\Documents and Settings\user\Local Settings\Temp\ST_0200003t
> .tmp , replace
invalid 'and' 
      --------------------------------------------- end akdensity0.ipolate4 ---
      qui means `fipol' [aw=`tt'] if `use'
      qui gen `double' `lambda' = sqrt(r(mean_g)/`fipol') if `use'
    -------------------------------------------------------- end akdensity0 ---
    qui gen double `tmph' = (`wwidth')*(`tmplambda')
    di as text "Step 2: Adaptive kernel density estimation"
    qui akdensity0 `ix' [`weight'`exp'] if `use' , at(`m') generate(`d') width(
> `tmph') `nbands' double
  ----------------------------------------------------------- end akdensity ---

My code runs as follows:

/*calculating the bandwith using Deaton's (1997) formula in turn based on Silverman's (1986) rule of thumb*/

local fname1 "c:\data80\hhcons80"
local fname2 "c:\data85\hhcons85"
local fname3 "c:\data90\hhcons90"
local fname4 "c:\data95\hhcons95"
	local l=1
	while `l'<=4 {
		drop _all
		use `fname`l'', clear
		keep x hhsize
		qui gen xx=x
		qui gen n=hhsize
		collapse (sum) n (sd) sdx=x (iqr) inqr=xx [aweight=n]
		qui gen q75= 0.75*inqr
		if sdx<q75{
		local h`l'=((2.34*sdx[1])/(n[1]^(1/5)))/2
		else {
		local h`l'=((2.34*q75)/(n[1]^(1/5)))/2
		dis "Half bandwith`l' = `h`l''"

/*combining consumption data and weights for all years in one file*/

		use `fname`l'', clear
		qui keep x hhsize
		qui gen n`l'=hhsize
		rename x x`l'
		qui keep x`l' n`l'
		if `l'<2{
		save akd, replace
		merge using akd
		drop _merge
		save akd, replace
		local l=`l'+1

/*calculating and graphing akdensities*/
 		use akd, clear
		akdensity x4 [aweight=n4], epanechnikov nogr width(`h4') gen(y4 fy4)
		akdensity x3 [aweight=n3], epanechnikov nogr width(`h3') gen(fy3) at(y4)
		akdensity x2 [aweight=n2], epanechnikov nogr width(`h2') gen(fy2) at(y4)
		akdensity x1 [aweight=n1], epanechnikov nogr width(`h1') gen(fy1) at(y4) 
		label var fy4 "1995"
		label var fy3 "1990"
		label var fy2 "1985"
		label var fy1 "1980"

		line fy4 fy3 fy2 fy1 y4 if y4<2000, sort ytitle(Density) xtitle(Per capita monthly consumption (Rs))title( "Distribution of Per Capita Consumption, 1980-96") leg(rows(1)) clpattern(shortdash_dot longdash_dot dash solid) 
		graph save akd.gph, replace


Would anybody have any ideas re: what's going wrong?
Thanks so much,
*   For searches and help try:

© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index