************************************************
*! xtvariog2 1.1 AN 27 Oct 2009 added panel support
*! based on variog2 1.0.0 NJC 4 March 2005
program xtvariog2, sort
version 8
syntax varlist(numeric min=3 max=3) [if] [in], Width(numlist max=1 >0) ///
[ Lags(real 0) LIST Generate(string) Id(varlist) noGRaph * ]
if "`id'"=="" {
xtset
loc id=r(panelvar)
}
tempvar ivar t
g long `t'=_n
egen long `ivar'=group(`id')
* check data and options
marksample touse
markout `touse' `ivar'
qui count if `touse'
if r(N) == 0 error 2000
local N = r(N)
if "`generate'" != "" {
confirm new variable `generate'
conf new variable `generate'_lag
}
if `lags' < 0 {
di as err "lags must be positive"
exit 198
}
tempvar ict
egen long `ict'=sum(`touse'), by(`ivar')
su `ict' if `ivar'<., mean
loc maxn=r(max)
local lags = cond(`lags',`lags',int(`maxn'/2))
local lags = min(`lags',`maxn'-5)
* calculation
tempvar lag gamma npairs LAG
tokenize "`varlist'"
args z x y
quietly {
replace `touse' = -`touse'
sort `touse' `ivar' `t'
by `touse' `ivar': gen `lag' = _n if _n<`lags' & `touse'==-1
label var `lag' "Lag (width `width')"
gen byte `LAG' = .
gen double `gamma' = 0
gen long `npairs' = 0
label var `gamma' "Semi-variance"
forval i = 1/`maxn' {
local I = `i' + 1
forval j = `I'/`N' {
by `touse' `ivar':replace ///
`LAG'=ceil(sqrt((`y'[`i']-`y'[`j'])^2+((`x'[`i']-`x'[`j'])^2))/`width')
by `touse' `ivar':replace `npairs'=`npairs'+1 if _n==`LAG'
by `touse' `ivar':replace `gamma'=`gamma'+(`z'[`i']-`z'[`j'])^2 if _n==`LAG'
}
}
replace `gamma' = `gamma' / (2 * `npairs')
}
* list if desired
if "`list'" == "list" {
char `lag'[varname] "Lag"
char `gamma'[varname] "Semi-variance"
char `npairs'[varname] "# of pairs"
l `lag' `gamma' `npairs' if !mi(`lag',`gamma'), ///
subvarname noobs abb(13) sepby(`ivar')
}
* graph
if "`graph'"=="" {
local zlab : variable label `z'
if `"`zlab'"' == "" local zlab "`z'"
loc o ti(`"Semi-variogram of `zlab'"') i(`ivar') t(`lag') leg(off) `options'
xtline `gamma' if `gamma'<., ov ysc(r(0 .)) yla(, ang(h)) `o'
}
* generate if desired
qui if "`generate'" != "" {
ren `gamma' `generate'
label var `generate' "Semi-variance of `zlab'"
ren `lag' `generate'_lag
}
end
************************************************
use http://www.ats.ucla.edu/stat/stata/faq/ozone
g obs=_n
expand 2
bys station: g year=_n
drop if year==2 & obs>14
set seed 1
replace av8top=av8top+rnormal() if year==2
xtvariog2 av8top lat lon, w(.1) i(year) list g(x)
On Tue, Oct 27, 2009 at 10:00 AM, Stas Kolenikov <[email protected]> wrote:
> The standard recommendation applies: take the code of -variog*- and
> modify it to suit your needs (saving a copy under a different name).
> You would have to change all the single observation concepts in
> -variog*- with the panel-level concepts (which is moderately
> difficult), and add another cycle in the computation of `zdsq' scalar
> (which is quite easy).
>
> On 10/26/09, Raphael Fraser <[email protected]> wrote:
>> Are there any programs available to estimate variogram for
>> longitudinal data? Note that -variog- or -variog2- does not apply to
>> longitudinal data.
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
