Hi Bugs users,
I have tried to find some SAS programs to diagnose the convergence of BUGS output in SAS. But I didn't find this kind of scripts. Then I decide to write one. I chose the Geweke and Yu & Mykland methods. When calculting the Geweke, I used the spectral density from PROC SPECTRA which depends on the choice of kernel. I am not familiar with time series. Then I used the quadratic kernel. However, I got quite different answers with R-CODA. I attached the scripts with this email. And I need your help on how to choose the weight in proc spectra. I really appreciate any help you can provide on this.
Thanks,
Zhiyong Zhang
Department of Psychology
University of Virginia
Email: [log in to unmask]
/*These codes have been written to perform convergence diagnosis in SAS.
Please send questions or comments to Zhiyong Zhang ([log in to unmask]).
Two statistics were calculated - Geweke (1992) and Yu and Mykland (1994) with
D statistics (Brook, 1996, cf. Brooks & Roberts, 1995).
Geweke falls in (-1.96, 1.96) and D falls in (.5 +/- Za/2*sqrt(1/4n)) mean convergence.
*/
options symbolgen;
/*Change the data and variable you want to test*/
%let datafile=post1;
%let varlist=mul;
/*Get the # of observations of the dataset specified*/
%let dsid = %sysfunc(open(&datafile));
%let nobs =%sysfunc(attrn(&dsid,NOBS));
%let rc = %sysfunc(close(&dsid));
/* Data for the first 10% */
%let leg10=&nobs;
%let leg10=&leg10*0.1;
%let leg10=%sysevalf(&leg10);
data first10;
set &datafile(firstobs=1 obs=&leg10);
run;
/* Data for the last 50% */
%let leg50=&nobs;
%let leg50=&leg50*0.5;
%let leg50=%sysevalf(&leg50+1);
data last50;
set &datafile (firstobs=&leg50);
run;
/*Calculate the mean for the first 10%*/
proc means data=first10 noprint;
var &varlist;
output out=m10 mean=avg;
run;
data _null_;
set m10;
call symput('avg10',avg);
run;
*%put &avg10;
/*Calculate the mean for the last 50%*/
proc means data=last50 noprint;
var &varlist;
output out=m50 mean=avg50;
run;
data _null_;
set m50;
call symput('avg50',avg50);
run;
*%put &avg50;
/*Calculate the spectral density*/
/*For first 10%*/
proc spectra data=first10 out=b10 p s;
weight qs;
var &varlist;
run;
data b10;
set b10 (firstobs=1 obs=1);
call symput('var10',s_01);
run;
*%put &var10;
/*For last 50%*/
proc spectra data=last50 out=b50 p s;
weight qs;
var &varlist;
run;
data b50;
set b50 (firstobs=1 obs=1);
call symput('var50',s_01);
run;
*%put &var50;
/*the function to calculate the Geweke statistics*/
%macro geweke(avg10, avg50, var10, var50, n1, n2);
(&avg10-&avg50)/sqrt(1/&n1*&var10+1/&n2*&var50);
%mend;
data temp;
geweke=%geweke(&avg10,&avg50,&var10,&var50,&leg10,&leg50);
output;
run;
proc print data=temp;
run;
proc means data=&datafile noprint;
var &varlist;
output out=avg mean=avg;
run;
data _null_;
set avg;
call symput('avg',avg);
run;
*%put &avg;
data YuMyk;
set &datafile;
mu=&avg;
retain st nmu 0;
st=&varlist+st;
nmu=mu+nmu;
sthat=st-nmu;
N=_N_;
run;
proc gplot data=YuMyk;
plot sthat*N;
run;
data temp;
set YuMyk;
stminus1=lag2(&varlist);
stself=lag1(&varlist);
stplus1=&varlist;
keep stminus1 stplus1 stself;
run;
data temp;
set test (firstobs=3);
d=0;
if (stself<stminus1 & stself<stplus1) or (stself>stminus1 & stself>stplus1) then d=1;
run;
proc means data=temp;
var d;
run;
-------------------------------------------------------------------
This list is for discussion of modelling issues and the BUGS software.
For help with crashes and error messages, first mail [log in to unmask]
To mail the BUGS list, mail to [log in to unmask]
Before mailing, please check the archive at www.jiscmail.ac.uk/lists/bugs.html
Please do not mail attachments to the list.
To leave the BUGS list, send LEAVE BUGS to [log in to unmask]
If this fails, mail [log in to unmask], NOT the whole list
|