/*  This program generates realizations from a Stochastic Dummy Seasonality model ala
    A.C. Harvey (1989).
    This model turns out to be mathematically equivalent to an
    SMA(1) model on seasonal span differenced data.  */

options pagesize=60 linesize=74 nodate;                                                                                                    
goptions device=win                                                                                                                        
gsfname=plot                                                                                                                               
rotate=landscape                                                                                                                           
gsfmode=append                                                                                                                             
target=hpljs3;                                                                                                                             
*border; 
ddata mc;
     gamma_1=0;
	 gamma_2=0;
	 gamma_3=0;
 do t = -100 to 1000;
     omega = rannor(2982);
     epsilon = rannor(226);
     gammat = -gamma_1 - gamma_2 - gamma_3 + 0.01*omega;
	 y = gammat + epsilon;
	 if t > 0 then output;
     gamma_3 = gamma_2;
	 gamma_2 = gamma_1;
	 gamma_1 = gammat;
	 end;
      

data mc;
  set mc;
  d4y = y - lag4(y);
  
run;

proc print data = mc;

	run; 

title 'ACF and PACF for the Stochastic Dummy Seasonal model';
proc arima data=mc;
  identify var = y;

  run;

title 'ACF and PACF for Seasonal Span Differenced Stochastic Dummy Seasonal data';
proc arima data=mc;
  identify var = d4y;

  run;

proc gplot data=mc;                                                                                                                       
  symbol v=dot c=black i=join h=.8;                                                                                                        
  title1 'Stochastic Dummy Seasonal Model';                                                                                    
  title2 'X=Time Y=Stochastic Dummy Seasonal data';                                                                                                          
  axis1 order=(0 to 100 by 10)                                                                                                             
  label=(f=duplex 'Obs');                                                                                                                  
  axis2 order=(-25 to 25 by 5)                                                                                                              
  label=(f=duplex 'Y Series');                                                                                                         
  plot y*t / haxis=axis1 vaxis=axis2;
                                                                                                                                          
  run;  


proc gplot data=mc;                                                                                                                       
  symbol v=dot c=black i=join h=.8;                                                                                                        
  title1 'Seasonal Span Differenced stochastic dummy seasonal data';                                                                                    
  title2 'X=Time Y=D4Y Series';                                                                                                          
  axis1 order=(0 to 100 by 10)                                                                                                             
  label=(f=duplex 'Obs');                                                                                                                  
  axis2 order=(-25 to 25 by 5)                                                                                                              
  label=(f=duplex 'DY Series');                                                                                                         
  plot d4y*t / haxis=axis1 vaxis=axis2;
                                                                                                                                          
  run;