/* This program generates non-seasonal Airline data for comparison
   in the Buys Ballot Plots. The non-seasonal Airline data was generated
   by Monte Carlo with the same intercept, slope, and standard error of
   regression as that of the original log(passengers) data.  The following
   Buys Ballot Graphs are compared: "Seasons By Year", "Seasonal Trends",
   and "Seasonal Averages".  The comparative patterns are quite distinctive.
   Also we use the Deterministic Trend / Deterministic and Seasonal Model
   with AR(4) and a quadratic trend term (to allow for generality) to
   do a joint F-test for the presence of seasonality.  As can be seen,
   the joint F-test is highly significant in the original Airline
   passenger data while, as expected, it is insignificant in the Monte Carlo
   Non-seasonal Airline data. Also as a supplementary test for the presence
   or absence of seasonality in these two data sets we use Proc ARIMA to
   examine the ACF of the undifferenced and first differenced series to see
   the difference in their behaviors.  */
 
data airline;
   format date monyy5.;
   input date:monyy5. pass @@;
   cards;
jan49 112 feb49 118 mar49 132 apr49 129 may49 121 jun49 135
jul49 148 aug49 148 sep49 136 oct49 119 nov49 104 dec49 118
jan50 115 feb50 126 mar50 141 apr50 135 may50 125 jun50 149
jul50 170 aug50 170 sep50 158 oct50 133 nov50 114 dec50 140
jan51 145 feb51 150 mar51 178 apr51 163 may51 172 jun51 178
jul51 199 aug51 199 sep51 184 oct51 162 nov51 146 dec51 166
jan52 171 feb52 180 mar52 193 apr52 181 may52 183 jun52 218
jul52 230 aug52 242 sep52 209 oct52 191 nov52 172 dec52 194
jan53 196 feb53 196 mar53 236 apr53 235 may53 229 jun53 243
jul53 264 aug53 272 sep53 237 oct53 211 nov53 180 dec53 201
jan54 204 feb54 188 mar54 235 apr54 227 may54 234 jun54 264
jul54 302 aug54 293 sep54 259 oct54 229 nov54 203 dec54 229
jan55 242 feb55 233 mar55 267 apr55 269 may55 270 jun55 315
jul55 364 aug55 347 sep55 312 oct55 274 nov55 237 dec55 278
jan56 284 feb56 277 mar56 317 apr56 313 may56 318 jun56 374
jul56 413 aug56 405 sep56 355 oct56 306 nov56 271 dec56 306
jan57 315 feb57 301 mar57 356 apr57 348 may57 355 jun57 422
jul57 465 aug57 467 sep57 404 oct57 347 nov57 305 dec57 336
jan58 340 feb58 318 mar58 362 apr58 348 may58 363 jun58 435
jul58 491 aug58 505 sep58 404 oct58 359 nov58 310 dec58 337
jan59 360 feb59 342 mar59 406 apr59 396 may59 420 jun59 472
jul59 548 aug59 559 sep59 463 oct59 407 nov59 362 dec59 405
jan60 417 feb60 391 mar60 419 apr60 461 may60 472 jun60 535
jul60 622 aug60 606 sep60 508 oct60 461 nov60 390 dec60 432
;

/* Extract year and month series */

data airline;
   set airline;
   label series_yr = "Series Year"
         series_mon = "Month";

   series_yr = year(date);
   series_mon = month(date);

run;

/* Create dummy variables etc. */

data airline;
  set airline;
  lpass = log(pass);
  t = _n_;
  t2 = t*t;
  d1 = (series_mon=1);
  d2 = (series_mon=2);
  d3 = (series_mon=3);
  d4 = (series_mon=4);
  d5 = (series_mon=5);
  d6 = (series_mon=6);
  d7 = (series_mon=7);
  d8 = (series_mon=8);
  d9 = (series_mon=9);
  d10 = (series_mon=10);
  d11 = (series_mon=11);
  d12 = (series_mon=12);

run;

/* Determine the trend and standard error of regression to use for 
   non-seasonal Airline data. */

ods html;
  ods graphics on;

proc reg data = airline;
  model lpass = t;

  run;

/* Use Monte Carlo to generate non-seasonal Airline data. */ 

data mc;

do t = 1 to 144;
  a = rannor(440285);

lpassns = 4.81 + 0.01*t + 0.14*a;

  output;

end;

keep lpassns;

/* Combine the Monte Carlo data with original data for comparisons to follow. */

data combine;
   merge airline mc;

run;

/*Creates a Unique line color for each of the 12 years or months whichever is 
  called for.  */

symbol1 color=black interpol=spline
        value=none;
symbol2 color=blue interpol=spline
        value=none;
symbol3 color=brown interpol=spline
        value=none;
symbol4 color=libg interpol=spline
        value=none;
symbol5 color=cyan interpol=spline
        value=none;
symbol6 color=gold interpol=spline
        value=none;
symbol7 color=gray interpol=spline
        value=none;
symbol8 color=green interpol=spline
        value=none;
symbol9 color=lilac interpol=spline
        value=none;
symbol10 color=steel interpol=spline
        value=none;
symbol11 color=magenta interpol=spline
        value=none;
symbol12 color=maroon interpol=spline
        value=none;

/* The major graphs follow after this.  */

title1 'Log of Airline Passengers Jan. 1949 - Dec. 1960';
title2 '(Log of thousands of passengers)';

axis1 label=('Year');
axis2 order=(4 to 7 by .25)
      label=(angle=90 'Log of Passengers');

proc gplot data=combine;
   plot lpass*date / haxis=axis1 vaxis=axis2;
   symbol1 i=join;
   format date year4.;

run;

title1 'Log of Non-Seasonal Airline Passengers Jan. 1949 - Dec. 1960';
title2 '(Log of thousands of passengers)';

axis1 label=('Year');
axis2 order=(4 to 7 by .25)
      label=(angle=90 'Log of Passengers');

proc gplot data=combine;
   plot lpassns*date / haxis=axis1 vaxis=axis2;
   symbol1 i=join;
   format date year4.;

run;

/*Sets the axis parameters*/

axis1 label=('Month') minor=none
      value=('Jan' 'Feb' 'Mar' 'Apr' 'May' 'Jun'
	  'Jul' 'Aug' 'Sep' 'Oct' 'Nov' 'Dec')
      order=(1 to 12 by 1)
      offset = (2)
      width = 3;
axis2 label=(angle=90 'Log Passengers')
      order=(4 to 7 by 0.2);

title1 'Log of Airline Passenger Volume 1949 to 1960';
title2 '(Log of Thousands of Passengers) ';
title2 'Buys Ballot "Seasons By Year" Plot';

proc gplot data=combine;
   plot lpass*series_mon=series_yr / haxis=axis1
        vaxis=axis2 vminor=1 legend=legend1;
run;

/*Sets the axis parameters*/

axis1 label=('Month') minor=none
      value=('Jan' 'Feb' 'Mar' 'Apr' 'May' 'Jun'
	  'Jul' 'Aug' 'Sep' 'Oct' 'Nov' 'Dec')
      order=(1 to 12 by 1)
      offset = (2)
      width = 3;
axis2 label=(angle=90 'Log Passengers')
      order=(4 to 7 by 0.2);

title1 'Log of Non-seasonal Airline Passenger Volume 1949 to 1960';
title2 '(Log of Thousands of Passengers) ';
title2 'Buys Ballot "Seasons By Year" Plot';

proc gplot data=combine;
   plot lpassns*series_mon=series_yr / haxis=axis1
        vaxis=axis2 vminor=1 legend=legend1;
run;

/*Sets the axis parameters*/

axis1 label=('Year') minor=none
      value=('1949' '1950' '1951' '1952' '1953' '1954'
	  '1955' '1956' '1957' '1958' '1959' '1960')
      order=(1949 to 1960 by 1)
      offset = (2)
      width = 3;
axis2 label=(angle=90 'Log of Passengers')
      order=(4 to 7 by 0.2);

/*Graph the data - Seasonal Trends */

legend1 position = (bottom center outside);

title1 'Log of Airline Passenger Volume 1949 to 1960';
title2 '(Log of Thousands of Passengers) ';
title2 'Buys Ballot "Seasonal Trends" Plot';

proc gplot data=combine;
   plot lpass*series_yr=series_mon / haxis=axis1
        vaxis=axis2 vminor=1;

run;

/*Sets the axis parameters*/

axis1 label=('Year') minor=none
      value=('1949' '1950' '1951' '1952' '1953' '1954'
	  '1955' '1956' '1957' '1958' '1959' '1960')
      order=(1949 to 1960 by 1)
      offset = (2)
      width = 3;
axis2 label=(angle=90 'Log of Passengers')
      order=(4 to 7 by 0.2);

/*Graph the data - Seasonal Trends*/

legend1 position = (bottom center outside);

title1 'Log of Non-Seasonal Airline Passenger Volume 1949 to 1960';
title2 '(Log of Thousands of Passengers) ';
title2 'Buys Ballot "Seasonal Trends" Plot';

proc gplot data=combine;
   plot lpassns*series_yr=series_mon / haxis=axis1
        vaxis=axis2 vminor=1;

run;

  /* Here we test for seasonality using the joint F-test on the seasonal dummies
     using the Deterministic Trend / Deterministic Season Model with AR(4) errors
     and a quadratic time trend for generality. */

  title 'Joint F-test of Seasonality on Original log of Airline data';
  
proc autoreg data = combine;
  model lpass = t t2 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12/ nlag = 4 
                     method=ml backstep slstay=0.05;
     test d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12;
run;

  title 'Joint F-test of Seasonality on Non-seasonal Airline data';

proc autoreg data = combine;
  model lpassns = t t2 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12/ nlag = 4 
                     method=ml backstep slstay=0.05;
     test d2, d3, d4, d5, d6, d7, d8, d9, d10, d11, d12;
run;

/* Here we generate the various ACFs to see how they help identify the presence
   or absence of seasonality. */
 
title 'ACF of lpass';
proc arima data = combine;
   identify var = lpass;
   run;

title 'ACF of lpass(1)';
proc arima data = combine;
   identify var = lpass(1);
   run;

title 'ACF of lpassns';
proc arima data = combine;
   identify var = lpassns;
   run;

title 'ACF of lpassns(1)';
proc arima data = combine;
   identify var = lpassns(1);
   run;

/* We need to use Proc Sort to prepare our data for the use of Proc Means to
   help us in the construction of the "Seasonal Averages" Buys Ballot Plots.  */
     
proc sort data = combine;
   by series_mon;
   run;

/* Means By Month Over Years for both lpass and lpassns data.  */

title 'Means By Month Over Years for both lpass and lpassns data';

proc means data = combine;
   var lpass lpassns;
   by series_mon;
   output out=results_means mean = meanlpass meanlpassns;
   run;

title1 'Log of Airline Passenger Volume 1949 to 1960';
title2 '(Log of Thousands of Passengers) ';
title2 'Buys Ballot "Seasonal Averages" Plot';

proc gplot data = results_means;
	plot meanlpass*series_mon;
	symbol interpol=none value=dot;

	run;

title1 'Log of Non-Seasonal Airline Passenger Volume 1949 to 1960';
title2 '(Log of Thousands of Passengers) ';
title2 'Buys Ballot "Seasonal Averages" Plot';

proc gplot data = results_means;
	plot meanlpassns*series_mon;
	symbol interpol=none value=dot;

	run;