/*  This program generates realizations from a stochastic level stochastic slope model ala
    A.C. Harvey's (1989) book on structural time series models.  
    This model turns out to be mathematically equivalent to an
    ARIMA(0,2,2) Box-Jenkins model.  */

options pagesize=60 linesize=74 nodate;                                                                                                    
goptions device=win                                                                                                                        
gsfname=plot                                                                                                                               
rotate=landscape                                                                                                                           
gsfmode=append                                                                                                                             
target=hpljs3;                                                                                                                             
*border; 

data mc;
  mu1 = 0;
  beta1 = 0;
  do t = 0 to 100;
     eta = rannor(494);
	 zeta = rannor(9937);
	 epsilon = rannor(647);
	 mu = mu1 + beta1 + eta;
	 beta = beta1 + zeta; 
     y = mu + 2.0*epsilon;
     if t > 0 then output;
     mu1 = mu;
	 beta1 = beta;
     end;

proc gplot data=mc;                                                                                                                       
  symbol v=dot c=black i=join h=.8;                                                                                                        
  title1 'Monte Carlo of the Stochastic Level Stochastic Slope model';                                                                                    
  title2 'X=Time Y=Stochastic Level Stochastic Slope data';                                                                                                          
  axis1 order=(0 to 100 by 10)                                                                                                             
  label=(f=duplex 'Obs');                                                                                                                  
  axis2 order=(-100 to 400 by 25)                                                                                                              
  label=(f=duplex 'Y');                                                                                                         
  plot y*t / haxis=axis1 vaxis=axis2;

data mc;
  set mc;
  dy = y - lag1(y);
  ddy = dy - lag(dy);

run;

title 'ACF and PACF for the local trend stochastic slope model';
proc arima data=mc;
  identify var = y;

  run;

title 'ACF and PACF for First Differenced Stochastic Level Stochastic Slope data';
proc arima data=mc;
  identify var = y(1);

  run;

title 'ACF and PACF for Second Differenced Stochastic Level Stochastic Slope data';
proc arima data=mc;
  identify var = y(2);

  run;

proc gplot data=mc;                                                                                                                       
  symbol v=dot c=black i=join h=.8;                                                                                                        
  title1 'First Differenced Stochastic Level Stochastic Slope data';                                                                                    
  title2 'X=Time Y=DDY Series';                                                                                                          
  axis1 order=(0 to 100 by 10)                                                                                                             
  label=(f=duplex 'Obs');                                                                                                                  
  axis2 order=(-40 to 40 by 5)                                                                                                              
  label=(f=duplex 'DY Series');                                                                                                         
  plot dy*t / haxis=axis1 vaxis=axis2;
                                                                                                                                          
  run;  



proc gplot data=mc;                                                                                                                       
  symbol v=dot c=black i=join h=.8;                                                                                                        
  title1 'Second Differenced Stochastic Level Stochastic Slope data';                                                                                    
  title2 'X=Time Y=DDY Series';                                                                                                          
  axis1 order=(0 to 100 by 10)                                                                                                             
  label=(f=duplex 'Obs');                                                                                                                  
  axis2 order=(-40 to 40 by 5)                                                                                                              
  label=(f=duplex 'DY Series');                                                                                                         
  plot ddy*t / haxis=axis1 vaxis=axis2;
                                                                                                                                          
  run;