***************************************************; /*Parameters to be filled:*/ /*Baseint - the estimated baseline intercept from the final model*/ /*Baseslope - the estimated baseline slope from the final model*/ /*Interint - the estimated level change comparing the intervention period to the baseline from the final model*/ /*Interslope - the estimated trend change comparing the intervention period to the baseline from the final model */ /*nperiod - the number of total time points*/ /*Sigmas - the estimated variance from the final model*/ /*Ar - the estimated autocorrelation terms from the final model */ /*Iter - the number of simulation iterations */ /*Intm - the value of the starting time point of the first policy*/ /*Modar - the maximum allowable autocorrelations will be used in each iteration of segmented regression desired by the investigators*/ /*Phaseup - the upper bound of the phase-in period for the policy*/ /*Phaselow - the lower bound of the phase-in period for the policy*/ /*Alpha- the significance level, usually it is 0.05*/ /*Fieldb- the list of all terms from the full model*/ /*Fields- the list of all terms except the terms of interest, policy terms*/ /*Numdif - the number of degree freedoms*/ /*Outlog - the output file including power*/ /*Batch - the number of simulations one wants to conduct simultaneously*/ ****************************************************; options ps=20000 nonotes nosyntaxcheck nosymbolgen; libname tmz 'c:\'; %macro sim_oneint_onegrp (baseint,baseslope,interint,interslope,nperiod,sigmas,ar,iter,intm,modar,phaseup,phaselow, alpha,fieldb,fields,numdif,outlog,batch); %macro numa;%do i=1 %to &iter;logd&i%end;%mend numa; data xbeta(drop=beta0-beta3); beta0=&baseint;beta1=&baseslope;beta2=&interint;beta3=&interslope; do trend=1 to &nperiod;lvchg=(trend>=&intm); if lvchg=0 then trchg=0;else trchg=trend+1-&intm; mu2=beta0+beta1*trend+beta2*lvchg+beta3*trchg; output; end; run; data xbetab(drop=i); set xbeta; do i=1 to &batch; simby=i; output; end; run; proc sort data=xbetab; by simby trend; run; %do k=1 %to &iter; proc iml symsize=309600000 worksize=3000000000; reset noname; use xbetab var {mu2}; read all into z2; start armasim(y,n,phi,theta,seed,sigma_2); p = ncol(phi)-1;q = ncol(theta)-1; y = normal(j(1,n+q,seed));y = sqrt(sigma_2)*y; /*---- Pure MA or white noise ----*/ if p=0 then y=product(theta,y)[,(q+1):(n+q)]; else do; /*---- Pure AR or ARMA ----*/ /*---- Get the autocovariance function ----*/ call armacov(gamma,cov,ma,phi,theta,p); if gamma[1]<0 then do; print 'ARMA parameters not stable.'; print 'Execution terminating.'; stop; end; /*---- Form covariance matrix ----*/ gamma=toeplitz(gamma); /*---- Generate covariance between initial y and ----*/ /*---- initial innovations ----*/ if q>0 then do; psi = ratio(phi,theta,q); psi = hankel( psi[,-( (-q):(-1) ) ] ); m = max(1,(q-p+1)); psi = psi[ -((-q):(-m)) ,]; if p>q then psi = j(p-q,q,0)//psi; gamma = ( gamma||psi ) // ( psi`||i(q) ); end; /*---- Use Cholesky root to get startup values ----*/ gamma = root(gamma); startup = y[,1:(p+q)]*gamma; e = y[,(p+q+1):(n+q)]; /*---- Generate MA part ----*/ if q>0 then do; e = startup[,(p+1):(p+q)]||e; e = product(theta,e)[,(q+1):(n+q-p)]; end; y = startup[,1:p]; phi1 = phi[, -(-(p+1):(-2)) ]`; /*---- Use difference equation to generate ----*/ /*---- remaining values ----*/ do ii=1 to n-p; y = y||( e[,ii]-y[,ii:(ii+p-1)]*phi1 ); end; end; y = y`; finish; /*---- ARMASIM ----*/ seed=%eval(&k*&iter); sample=%eval(&nperiod*&batch+200); sam=sample+1-%eval(&nperiod*&batch); run armasim(y,sample,{&ar},{1},seed,&sigmas); ya=y[sam:sample]; yf=ya+z2; create ts from yf; append from yf; quit; filename routed 'c:\zzzz.txt'; proc printto print=routed new; run; data dd;merge xbetab ts; if &phaseup ne . then do;if &phaseup>=trend>=&phaselow then col1=.;end; else do;end; run; %if &modar=0 %then %do; proc autoreg data=dd; ods select FitSummary; ods output FitSummary=work.logk1; model col1=&fieldb/method=ml loglikl; by simby; run; data t1;set logk1(keep=nvalue1 Label1 simby rename=(nvalue1=logk1) where=(Label1='Log Likelihood')); run; proc autoreg data=dd; ods select FitSummary; ods output FitSummary=work.logk2; model col1=&fields/method=ml loglikl; by simby; run; data t2;set logk2(keep=nvalue1 Label1 simby rename=(nvalue1=logk2) where=(Label1='Log Likelihood')); run; data logd(drop=label1);merge t1 t2;by simby;ldiff=-2*(logk2-logk1);run; %end; %else %do; proc autoreg data=dd; ods select FitSummary; ods select ConvergenceStatus; ods output FitSummary=work.logk1; ods output ConvergenceStatus=work.conv1; model col1=&fieldb/method=ml nlag=&modar loglikl; by simby; run; data t1;set logk1(keep=nvalue1 Label1 simby rename=(nvalue1=logk1) where=(Label1='Log Likelihood')); run; data t1;set t1;by simby;if last.simby;run; proc autoreg data=dd; ods select FitSummary; ods select ConvergenceStatus; ods output FitSummary=work.logk2; ods output ConvergenceStatus=work.conv2; model col1=&fields/method=ml nlag=&modar loglikl; by simby; run; data t2;set logk2(keep=nvalue1 Label1 simby rename=(nvalue1=logk2) where=(Label1='Log Likelihood')); run; data t2;set t2;by simby;if last.simby;run; data logd;merge t1 t2;by simby;ldiff=-2*(logk2-logk1);run; data conv;merge conv1(keep=status simby) conv2(keep=status simby rename=(status=sta)); by simby; maxc=max(status,sta);run; data logd(drop=label1);merge logd conv(keep=maxc simby);by simby;run; %end; data logd;set logd;calpha=cinv(1-&alpha,&numdif);if maxc=0 then power=(ldiff>calpha);else power=.;run; proc printto;run; proc append base=logall data=logd;run; %end; data logall;set logall;jj=1;run; proc sql; create table outlog as select jj,sum(power=0) as p0,sum(power=1) as p1 from logall(firstobs=2) group by jj order by jj; quit; data tmz.&outlog;set outlog;power=p1/(p0+p1);run; %mend sim_oneint_onegrp; data logall;calpha=0;logk1=0;logk2=0;ldiff=0;power=0;maxc=0;simby=0;run; %sim_oneint_onegrp(0.6,0.0,0.00002,0.0,180,0.0000000004,%str(1 -0.5),10,73,1,0,0,0.05,%str(lvchg),%str(),1,pw1,10);