disp(['Running MATLAB script s133egstrikes.m']) ;
%
%  Makes Plot about "stationarity" for strikes data
%  show real data, together with simulated data
%  (created by old Gauss program 133eg1.prg)



%  read in Brockwell Davis data set
%
disp('   Reading Strikes data') ;
infstr = 'strikes.tsm' ;
fid = fopen(infstr,'r') ;
  strikedata = fscanf(fid,'%g') ;
fclose(fid) ;



%  set parameters for strikes data
%
%size(strikedata) 
n = length(strikedata) ;
sdrange = max(strikedata) - min(strikedata) ;
sdbottom = min(strikedata) - 0.10 * sdrange ;
sdtop = max(strikedata) + 0.10 * sdrange ;



%  Plot strikes
%
subplot(2,2,1) ;
  plot((1951:1980)',strikedata,'-o') ;
    title('Strikes Data') ;
    axis([1951,1980,sdbottom,sdtop]) ;
    xlabel('year') ;
    ylabel('# of strikes') ;




%  read in simulated data
%
disp('   Reading Simulated data') ;
infstr = '133egstr.out' ;
fid = fopen(infstr,'r') ;
  [simdata, cnt] = fscanf(fid,'%g') ;
fclose(fid) ;

simdata = reshape(simdata,4,cnt/4)' ;
          %   transpose, because of way input goes



vt = simdata(:,1) ;
          %  vector of "times"
vx = simdata(:,2) ;
          %  vector of "x_s"
vmu = simdata(:,3) ;
          %  vector of "meanss"
vssig = simdata(:,4) ;
          %  vector of "scaled sds"



%  set parameters
%
simleft = 0 ;
simright = 500 ;
simmin = min(vx) ;
simmax = max(vx) ;
simrange = simmax - simmin ;
simmin = simmin - 0.10 * simrange ;
simmax = simmax + 0.10 * simrange ;



%  plot big process
%
subplot(2,2,4) ;
  plot(vt,vx,'-') ;
    title('Stationary X_t = MU_t + SIG_t * N(0,1)') ;
    axis([simleft,simright,simmin,simmax]) ;
    xlabel('t') ;
    ylabel('X_t') ;



%  Set parameters for chopped process
%
chopleft = 173 ;
chopright = 225 ;
chopind = find(chopleft <= vt  &  vt <= chopright) ;
chopt = vt(chopind) ;
chopx = vx(chopind) ;
chopmin = min(chopx) ;
chopmax = max(chopx) ;
choprange = chopmax - chopmin ;
chopmin = chopmin - 0.10 * choprange ;
chopmax = chopmax + 0.10 * choprange ;



%  add lines to big process plot
%
hold on ;
  plot([chopleft; chopleft],[simmin,simmax]) ;
  plot([chopright; chopright],[simmin,simmax]) ;
hold off ;



%  plot chopped process
%
subplot(2,2,2) ;
  plot(chopt,chopx,'-o') ;
    title('Simulated, Stationary X_t') ;
    axis([chopleft,chopright,chopmin,chopmax]) ;
    xlabel('t') ;
    ylabel('X_t') ;



%  plot driving processes
%
drivmin = 0 ;
drivmax = max(vx) ;
drivrange = drivmax - drivmin ;
drivmax = drivmax + 0.10 * drivrange ;
subplot(2,2,3) ;
  plot(vt,vmu,'-',vt,vssig,'--') ;
    title('Stationary Drivers, MU_t & SIG_t') ;
    axis([simleft,simright,drivmin,drivmax]) ;
    xlabel('t') ;
    ylabel('MU_t, 5*SIG_t') ;








  orient landscape ;
print \matlab\steve\ps\s133egstrikes.ps -dps ;



