/* AR(1) model */

proc import datafile="G:\Teaching\ECN410\ECN410.Beamers\Lecture11_ARMA\pepsi.price.06to16.xlsx" dbms=xlsx 
	out=work.pepsi replace; 
getnames=yes;
run;

data pepsi; set pepsi (firstobs=2);
run;

/* sort the data by date */
proc sort data=pepsi out=pepsi_sorted; by date;

data pepsi_sorted; set pepsi_sorted;
lagprice=lag(close);  /* get the lagged value of the closing price*/

proc print data=pepsi_sorted (obs=20);
run;

/* AR(1) model */
proc reg data=pepsi_sorted;
model close=lagprice;
run;
quit;



proc forecast data=pepsi_sorted 
out=es_pepsi outfull outest=estimates_pepsi method=expo  
nstart=1 trend=1 weight=.5 lead=1;
id date;
var close;
run;

proc sql; create table es_pepsi1 as select a.*, b.close as esf_close
	FROM pepsi as a LEFT JOIN es_pepsi as b
	ON a.date=b.date
	WHERE b._TYPE_="FORECAST"
	;
	QUIT;

data es_pepsi1; set es_pepsi1;
label esf_close = "forecast alpha=.5";
run;

symbol1 interpol=join;
symbol2 interpol=join;
legend1 label=none frame;

axis2 label=("Y");

proc gplot data=es_pepsi1;
plot (close esf_close)*date / overlay legend=legend1 haxis=axis1 vaxis=axis2;
run;
quit;





proc import datafile="G:\Teaching\ECN410\ECN410.Beamers\Lecture11_ARMA\XYZSALES11.xlsx" dbms=xlsx 
	out=work.xyzsales replace; 
getnames=yes;
run;

%let nperiods=3; /*choose the number of periods for your moving average*/
%let firstobs=%eval(&nperiods );

data xyzsales; set xyzsales;
time_period=_n_;


/* create moving average */
proc expand data=xyzsales out=movingaverage;
convert sales = sales_ma /transformout=(movave &nperiods);
run;

data movingaverage1; set movingaverage (firstobs=&firstobs);

run;

proc sql; create table movingaverage2 as select a.*, b.sales_ma as forecast
	FROM movingaverage as a LEFT JOIN movingaverage1 as b
	ON a.time_period=b.time_period+1
	;
	quit;

proc print data=movingaverage2;
var time_period sales /*sales_ma*/ forecast;
run;

data movingaverage2; set movingaverage2;
label forecast = "MA forecast";
SD=(sales-forecast)**2;
AD=abs(sales-forecast);
APE=abs(AD/sales);

proc means data=movingaverage2;
run;


/* create plot */

symbol1 interpol=join;
symbol2 interpol=join;
legend1 label=none frame;

axis2 label=("Y");

proc gplot data=movingaverage2;
plot (sales forecast)*time_period / overlay legend=legend1 haxis=axis1 vaxis=axis2;
run;
quit;



/* exponential smoothing */
proc forecast data=movingaverage 
out=exponential_smooth outfull /*outresid*/ 
outest=estimates
method=expo /* winters, expo*/ 
nstart=1
trend=1 /*2 for double exponenetial smoothing*/
weight=.3 
lead=1
;
id time_period;
var sales;

run;

proc print data=exponential_smooth;
where _TYPE_="FORECAST";
run;

proc print data=estimates;
run;

proc sql; create table exponential_smooth1 as select a.*, b.sales as esf_sales
	FROM movingaverage2 as a LEFT JOIN exponential_smooth as b
	ON a.time_period=b.time_period
	WHERE b._TYPE_="FORECAST"
	;
	QUIT;

data exponential_smooth1; set exponential_smooth1;
label esf_sales = "forecast alpha=.3";
run;

symbol1 interpol=join;
symbol2 interpol=join;
symbol3 interpol=join;
legend1 label=none frame;

axis2 label=("Y");

proc gplot data=exponential_smooth1;
plot (sales forecast esf_sales)*time_period / overlay legend=legend1 haxis=axis1 vaxis=axis2;
run;
quit;


/* double exponential smoothing for case where trend is present */

proc import datafile="G:\Teaching\ECN410\ECN410.Beamers\Lecture11_ARMA\NEWCON11.xlsx" dbms=xlsx 
	out=work.newcon replace; 
getnames=yes;
run;

data newcon_bench; set newcon;
trend=year-1990;

proc reg data=newcon_bench;
model newcon=trend;
run;
quit;

proc forecast data=newcon 
out=exponential_smooth2 outfull outest=estimates2
method=winters /* winters, expo*/ 
/*nstart=1*/
trend=2 /*2 for double exponenetial smoothing*/
weight=(.78 .23) 
lead=1
;
id year;
var newcon;

run;

proc print data=exponential_smooth2; 
run;

proc print data=estimates2;
run;



proc sql; create table exponential_smooth2a as select a.*, b.newcon as desf_newcon
	FROM newcon as a LEFT JOIN exponential_smooth2 as b
	ON a.year=b.year
	WHERE b._TYPE_="FORECAST"
	;
	QUIT;

data exponential_smooth2a; set exponential_smooth2a;
label desf_newcon = "forecast newcon";
run;

symbol1 interpol=join;
symbol2 interpol=join;
legend1 label=none frame;

axis2 label=("Y");

proc gplot data=exponential_smooth2a;
plot (newcon desf_newcon)*year / overlay legend=legend1 haxis=axis1 vaxis=axis2;
run;
quit;