#include <stdio.h>
#include <math.h>
#include <ctype.h>

#define	MAXOBS	150		/* Max number of observations */
#define	FALSE	0
#define TRUE	!FALSE
#define OK 	TRUE
#define ESC	27
#define FOREVER	TRUE
#define LINESPS 24		/* Number of lines per screen */

enum methods {MOVAVG, EXPSMO, DSMOOTH, TREND, MEAN, DECOMP, NONE};

int	begin;			/* Beginning period to print */
char	buf[80];		/* Scratch area */
char	*f1 = "%11.3lf";        /* Formats for reports */
char	*f2 = "%11s";
double	ferr[MAXOBS+1];		/* Forecast error for period t */
int	finish;			/* Last period to print is TSD */
double	fore[MAXOBS+1]; 	/* Forecast for "next" period */
double	infin = 1e100;		/* Infinity, for TSD */
int	lines;			/* Number of displayed lines on screen */
int	numv;			/* Number of observations */
int	pauseon = TRUE;		/* Is full-screen pause used? */
int	periods;		/* Number of observations/year */
char	seas_name[13][10];	/* Name of "season" */
int	start;			/* Number of obs in MSE calc */
double	y[MAXOBS+1];		/* Array of observations */

double	ma[MAXOBS+1];		/* Moving averages */
double	cma[MAXOBS+1];		/* Centered moving averages */
double	trend[MAXOBS+1];	/* trend for TSD */
double  cyc[MAXOBS+1];
double	rma[MAXOBS+1];
double	seas[13];
double	irreg[MAXOBS+1];
double	seas_hi[13];
double	seas_lo[13];
int	num_seas[13];

#define iseas(k) = (((k+periods-1) % periods)+1)

main(argc,argv)
int argc; char *argv[];
	{
	enum methods	choice;		/* Menu choice */
	enum methods	menu();		/* Menu choice selector */
	int	ok;

	greet();
	if (init_read(argc,argv) == !OK) return;
	while ((choice = menu()) != NONE)
		{
		switch(choice)
			{
			case MOVAVG:	ok = movavg(); break;
			case EXPSMO:	ok = expsmo(); break;
			case DSMOOTH:	ok = dsmooth(); break;
			case TREND:	ok = regress();  break;
			case MEAN: 	ok = mean();	  break;
			case DECOMP: 	ok = decomp();  break;
			}
		if (ok) summary();
		pause();
		}
	}

greet()
	{
	clrscr(); 
	cprintf("*** SMOOTHIE: INTERACTIVE TIME-SERIES FORECASTING ***\n");
	}

init_read(argc,argv)
int argc; char *argv[];
	{
	char	*filename = "FCST.DAT        ";
	FILE	*fn, *fopen();
	register int i;

	cprintf("by Richard S. Barr and Jim Collins\n");

	if (argc == 1) { help(); return !OK; }

	strcpy(filename,argv[1]);
	printf("\nProcessing data file  %s\n",filename);
	if ((fn = fopen(filename,"r")) == NULL)
		{
		printf("\n>> Could not find file.\n");
		return !OK;
		}
	for(i=1; i<MAXOBS; ++i)
		{
		if (fscanf(fn,"%lf",&y[i]) != 1) break;
		else	++numv;
		}
	printf("\n%d observations were read\n",numv);

	printf("\nHow many observations to use in the mean error calculation\n");
	printf("  (enter 0 for the maximum possible)? ");
	scanf("%d",&start);

	if(start<=0 || start>numv) start=numv;
	else start = numv-start+1;

	return OK;
	}

enum methods menu()
	{
	char	ch;

	greet();
	printf("\nThe following forecasting models are available:\n\n");
	cprintf("+----------------------------------+");
	cprintf("|  Moving average                  |");
	cprintf("|  Exponential smoothing (single)  |");
	cprintf("|  Smoothing with trend (double)   |");
	cprintf("|  Trend                           |");
	cprintf("|  Arithmetic mean                 |");
	cprintf("|  Decomposition                   |");
	cprintf("+----------------------------------+");

	printf("\nSelect model by typing the first letter, or Q to quit: ");

	while (TRUE)
		{
		ch=getchoice();
		switch(ch)
			{
			case 'M':	clrscr(); return MOVAVG;
			case 'E':	clrscr(); return EXPSMO;
			case 'S':	clrscr(); return DSMOOTH;
			case 'T':	clrscr(); return TREND;
			case 'A':	clrscr(); return MEAN;
			case 'D':	clrscr(); return DECOMP;
			case 'Q':	clrscr(); return NONE;
			default:	printf("Unknown choice %c",ch);
			}
		}
	}

cprintf(buf)
char *buf;
	{
	static char *blanks = "                                    ";
	register int indent;

	indent = (70-strlen(buf))/2;
	blanks[indent] = '\0';
	printf("%s%s\n",blanks,buf);
	blanks[indent] = ' ';
	}
movavg()
	{
	register int	i,k;
	int	periods;

	printf("             *** MOVING AVERAGE METHOD ***\n\n");

	for(periods=0;FOREVER;)
		{
		printf("Number of periods for your moving average: ");
		gets(buf);
		sscanf(buf,"%d",&periods);
		if (periods>0 && periods<=numv) break;
		printf(">> Periods must be > 0 and < %d\n",numv+1);
		}

	for(i=periods; i<=numv; ++i)
		{
		for (fore[i]=0.0, k=i-periods+1; k<=i; ++k) fore[i] += y[k];
		fore[i] /= ((double) periods);
		}

	begin	= periods;
	finish	= numv;
	lines	= 3;
	return TRUE;
	}

expsmo()
	{
	double alpha;		/* smoothing constant */
	register int i;

	printf("            *** EXPONENTIAL SMOOTHING (SINGLE) *** \n\n");
	while(FOREVER)
		{
		printf("What smoothing constant do you wish to use (0 < alpha < 1)? ");
		gets(buf);
		sscanf(buf,"%lf",&alpha);
		if (alpha>0.0 && alpha < 1.0) break;
		}

	fore[1]	= y[1];
	for (i=2; i<=numv; ++i) fore[i] = alpha*y[i] + (1.0-alpha)*fore[i-1];
	begin	= 1;
	finish	= numv;
	lines	= 3;
	return TRUE;
	}

dsmooth()
	{
	double alpha, beta;		/* smoothing constants */
	double snew, sold, tnew, told;
	register int i;

	printf("        *** EXPONENTIAL SMOOTHING WITH TREND (DOUBLE) *** \n\n");
	while(FOREVER)
		{
		printf("What DATA smoothing constant do you wish to use (0 < alpha < 1)? ");
		gets(buf);
		sscanf(buf,"%lf",&alpha);
		if (alpha>0.0 && alpha < 1.0) break;
		}
	while(FOREVER)
		{
		printf("What TREND smoothing constant do you wish to use (0 < beta < 1)? ");
		gets(buf);
		sscanf(buf,"%lf",&beta);
		if (beta>0.0 && beta< 1.0) break;
		}

	sold 	=  fore[1] = y[1];
	told	= 0.0;
	for (i=2; i<=numv; ++i)
		{
		snew	= alpha*y[i] + (1.0-alpha)*(sold+told);
		tnew	= beta*(snew-sold) + (1.0-beta)*told;
		fore[i] = snew + tnew;
		sold	= snew;
		told	= tnew;
		}

	printf("\nForecast using %lf + %lf * (number of periods beyond %d)\n",
		(float) sold, (float) told, numv);
	begin	= 1;
	finish	= numv;
	lines	= 6;
	return TRUE;
	}

mean()
	{
	double sumy;
	register int i;

	printf("             *** ARITHMETIC MEAN AS A FORECAST ***\n\n");

	for (sumy=0.0,i=1; i<=numv; ++i)
		{
		sumy += y[i];
		fore[i]	= sumy/i;
		}
	printf("The overall mean of Y is %lf\n",sumy/numv);
	begin	= 1;
	finish	= numv;
	lines	= 3;
	return TRUE;
	}

regress()
	{
	register int i;
	double sumx=0.0, sumy=0.0, sumxx=0.0, sumyy=0.0, sumxy=0.0;
	double	meanx,meany,a,b,rsq;
	printf("             *** LEAST-SQUARES TREND LINE ***\n\n");

	for(i=1; i<=numv; ++i)
		{
		sumx 	+= i;
		sumy 	+= y[i];
		sumxx	+= i*i;
		sumyy	+= y[i]*y[i];
		sumxy	+= i*y[i];
		meanx	= sumx/i;
		meany	= sumy/i;
		if (i>1)
			{
			b	= (sumxy-i*meanx*meany)/(sumxx-i*meanx*meanx);
			a	= meany-b*meanx;
			fore[i]	= a+b*(i+1);
			}
		else	fore[i]	= y[i];
		}
	rsq	= (a*sumy+b*sumxy-numv*meany*meany)/(sumyy-numv*meany*meany);
	printf("Trend equation:  T = %f + %f * t\n",(float)a,(float)b);
	printf("Sample coefficient of determination (R-squared) = %lf\n\n",rsq);
	begin	= 1;
	finish	= numv;
	lines	= 5;
	return TRUE;
	}

summary()
	{
	double	ferr, mse, mse2, sumerr=0.0, sum2=0.0;
	double	me=0.0, me2=0.0, mad=0.0, mad2=0.0, mape=0.0, mape2=0.0;
	int	numerr=0, num2=0;
	register int i;

	printf("\n                *** SUMMARY REPORT ***\n\n");
	printf("Period         Actual         Forecast          Error\n");
	printf("Number        Observn         for Next          (Forecast-\n");
	printf("                              Period            Actual)\n");
	printf("-----         -------         --------          --------\n");
	lines	+= 6;

	for(i=1; i<=numv; ++i)
		{
		printf(" %3d       %10.3lf",i,y[i]);
		if (i >= begin)
			{
			printf("      %10.3lf",fore[i]);
			if (i>begin)
				{
				ferr	= y[i]-fore[i-1];
				sumerr	+= ferr*ferr;
				mad	+= fabs(ferr);
				me	+= ferr;
				mape	+= fabs( 100.0*ferr/y[i] );
				++numerr;
				if (i >= start)
					{
					sum2 += ferr*ferr;
					mad2 += fabs(ferr);
					me2  += ferr;
					mape2 += fabs( 100.0*ferr/y[i] );
					++num2;
					}
				printf("        %10.3lf",ferr);
				}
			}
		printf("\n");
		if (++lines % (LINESPS-2) == 0) pause();
		}

	mse	= sumerr/numerr;
	me	= me/numerr;
	mad	= mad/numerr;
	mse2	= sum2/num2;
	me2	= me2/num2;
	mad2	= mad2/num2;
	mape	/= numerr;
	mape2	/= num2;

	if (LINESPS-lines <10) pause();

	printf("\n              *** FORECAST ERROR MEASURES ***\n\n");
	printf("                        \t%15s%15s\n","Specified","Maximum");
	printf("Number of periods       \t%15d%15d\n",num2,numerr);
	printf("Mean squared error (MSE)\t%15.3lf%15.3lf\n",mse2,mse);
	printf("Mean absolute error (MAD)\t%15.3lf%15.3lf\n",mad2,mad);
	printf("Mean absolute pct error (MAPE)\t%15.3lf%15.3lf\n",mape2,mape);
	printf("Mean error/Forecast bias\t%15.3lf%15.3lf\n",me2,me);

	}

decomp()
	{
	double back_ma, sumx=0.0, sumy=0.0, sumxx=0.0, sumyy=0.0, sumxy=0.0;
	double a, b, adj, rsq, sumirr=0.0, infin = 1e100, meanx, meany;
	int	num_tren,num_mod;
	register int i,j,k;

	printf("        *** TIME SERIES DECOMPOSITION ***\n\n");
	printf("Is the data monthly or quarterly (M/Q)? ");
	for(periods=0; periods == 0;)
		{
		switch (toupper(getchar()))
			{
			case 'M': printf("Monthly\n");
					periods = 12;
					break;
			case 'Q':	printf("Quarterly\n");
					periods = 4;
					break;
			}
		}


	set_seas_names();
			/* Set moving averages */
	back_ma	= periods/2 -1;
	for (i=periods; i<=numv; ++i)
		{
		k	= i-back_ma;
		for(ma[k]=0.0, j=(i-periods+1); j<=i; ++j) ma[k] += y[j];
		ma[k]	/= periods;
		}
			/* Centered moving averages */
	begin	= periods-back_ma;
	finish	= numv-back_ma-1;
	num_tren= finish-begin+1;
	if (num_tren < periods)
		{
		printf("\n>> There are insufficient data points to perform decomposition\n");
		pause();
		return FALSE;
		}
	for (i=begin; i<=finish; ++i)
		{
		cma[i]	= (ma[i]+ma[i+1])/2.0; 
		sumy	+= cma[i];
		sumx	+= i;
		sumxx	+= i*i;
		sumxy	+= i*cma[i];
		sumyy	+= cma[i]*cma[i];
		}
	meanx	= sumx/num_tren;
	meany	= sumy/num_tren;
	b	= (sumxy-num_tren*meanx*meany)/(sumxx-num_tren*meanx*meanx);
	a	= meany-b*meanx;
	rsq	= (a*sumy+b*sumxy-num_tren*meany*meany)/(sumyy-num_tren*meany*meany);
    printf("sumxy=%lf  meanx=%lf  meany=%lf   sumxx=%lf  num_tren=%d\n",sumxy,meanx,meany,sumxx,num_tren);

	printf("Trend equation:  T = %f + %f * Period\n",(float)a,(float)b);
	printf("Sample coefficient of determination (R-squared) = %lf\n",rsq);

			/* Trend and cyclical */
	for (i=begin; i<=finish; ++i)
		{
		trend[i]	= a+b*i;
		cyc[i]		= cma[i]/trend[i];
		}
			/* Ratio to MA*/
	for (i=1; i <=periods; ++i)
		{
		seas[i]		= 0.0;
		seas_hi[i] 	= -infin;
		seas_lo[i]	= infin;
		}
	for (i=begin; i<=finish; ++i)
		{
		rma[i]	= y[i]/cma[i];
		k	= ((i+periods-1) % periods)+1;
		seas[k]	+= rma[i];
		++num_seas[k];
		if (rma[i]>seas_hi[k]) seas_hi[k] = rma[i];
		if (rma[i]<seas_lo[k]) seas_lo[k] = rma[i];
		}
			/* Compute seasonal indices */
	for (i=1, num_mod=0, adj=0.0; i<=periods; ++i)
		{
		if (num_seas[i] >= 5)
			{
			seas[i]		-= (seas_hi[i]+seas_lo[i]);
			num_seas[i]	-= 2;
			++num_mod;
			}
		seas[i]	/= num_seas[i];
		adj	+= seas[i];
		}
			/* Adjust seasonal indices */
	printf("\nSeason       Seasonal Index\n");
	printf(  "------       --------------\n");
	adj = periods/adj;
	for(i=1; i<=periods; ++i)
		{
		seas[i]	*= adj;
		printf("%-10s       %8.4lf\n",&seas_name[i][0],seas[i]);
		}
	if (num_mod)
	  printf("\nModified means used in %d seasonal index calculations\n",num_mod);
	else
	  printf("\nUnmodified means used in seasonal index calculations\n");

	for (i=begin, sumirr=0.0; i<=finish; ++i)
		{
		irreg[i] = rma[i]/seas[((i+periods-1) % periods)+1];
		sumirr	+= irreg[i];
		}
	printf("Mean irregular index = %f\n",(float) sumirr/num_tren);

	printf("\nDo you want a detailed report ");
	if(yes()) dsummary(y,trend,seas,cyc);

	tsd_fcst(a,b,seas);
	return FALSE;

	}

tsd_fcst(a,b,seas)
double a,b,seas[];
	{
	int pd;
	double c,t,s,f;

	printf("\n            *** INTERACTIVE FORECASTING ***\n\n");

	while(FOREVER)
		{
		printf("Period to forecast or 0? ");
		scanf("%d",&pd);
		if (!pd) break;
		printf("Cyclical index to use or 1 to ignore? ");
		scanf("%lf",&c);
		if (c == 0.0) c = 1.0;
		t	= a+b*pd;
		s	= seas[((pd+periods-1) % periods)+1];
		f	= t * s * c;
		printf("   Forecast =       Trend * Seasonal Index * Cyclical Index\n");
		printf("%11.3lf = %11.3lf *  %11.3lf   * %11.3lf\n\n",f,t,s,c);
		}
	}

dsummary(y,trend,seas,cyc)
double y[],trend[],seas[],cyc[];
	{
	char	*dash=
	"------------------------------------------------------------------------";
	double	residual,mse;
	int	num_tren;
	register int j,k;

	printf("\n                 D E C O M P O S I T I O N    S U M M A R Y\n\n");
	printf("%s","Period");
	printf(f2,"Actual");
	printf(f2,"Trend");
	printf(f2,"Seasonal");
	printf(f2,"Cyclical");
	printf(f2,"Forecast");
	printf(f2,"Error");
	printf("\n");
	num_tren	= finish-begin+1;
	lines		= 4;
	for (j=begin, mse=0.0; j<=finish; ++j)
		{
		k	= (((j+periods-1) % periods)+1);
		fore[j]	= trend[j]*seas[k];
		residual	= y[j]-fore[j];
		mse	+= residual*residual;
		printf("  %4d",j);
		printf(f1,y[j]);
		printf(f1,trend[j]);
		printf(f1,seas[k]);
		printf(f1,cyc[j]);
		printf(f1,fore[j]);
		printf(f1,residual);
		printf("\n");
		if (!(j%periods)) { printf("%s\n",dash); ++lines;}
		if (((++lines) % (LINESPS-2)) == 0) pause();
		}
	printf("\nForecast based on trend and seasonal only\n");
	printf("Mean squared error for periods %d through %d = %f\n",
		begin,finish,(float)(mse/num_tren));
	pause();
	}


set_seas_names()
	{
	if (periods == 4)
		{
		register int i;
		for (i=1; i<=4; ++i)
			sprintf(&seas_name[i][0],"Quarter %d",i);
		return;
		}
	else
		{
		strcpy(&seas_name[1][0],"January");
		strcpy(&seas_name[2][0],"February");
		strcpy(&seas_name[3][0],"March");
		strcpy(&seas_name[4][0],"April");
		strcpy(&seas_name[5][0],"May");
		strcpy(&seas_name[6][0],"June");
		strcpy(&seas_name[7][0],"July");
		strcpy(&seas_name[8][0],"August");
		strcpy(&seas_name[9][0],"September");
		strcpy(&seas_name[10][0],"October");
		strcpy(&seas_name[11][0],"November");
		strcpy(&seas_name[12][0],"December");
		}
	}

yes()
	{
	static char ch;
	printf("(Y/N)? ");
	while (FOREVER)
		{
		gets(buf);
		if ((ch=buf[0]) == '\0') continue;
		if (islower(ch)) ch = toupper(ch);
		switch(ch)
			{
			case 'Y':	return TRUE;
			case 'N':	return FALSE;
			}
		}
	}

help()
	{
	char	*msg = "\
Program usage: 		smoothie filename\n\
\nwhere 'filename' is a file containing data for a single time series.\n\n\
*  The data are given in chronological order, with individual values separated\n\
   by blanks or placed on different lines.\n\
*  The program operates interactively, allowing the user to apply different\n\
   forecasting models to the same dataset.  Models available are:\n\n\
	- Moving averages\n\
	- Exponential smoothing, single and double\n\
	- Trend\n\
	- Arithmetic mean\n\
	- Full decomposition into trend, seasonal, and cyclical factors\n";

	printf("%s",msg);
	}
clrscr()
	{
	printf("");
	}

getchoice()
	{
	static char ch;
	do { gets(buf);}
	while (buf[0] == '\0');
	ch = buf[0];
	if (islower(ch)) ch = toupper(ch);
	return ch;
	}

pause()
	{
	static int longmessage = TRUE;
	if (!pauseon) return;
	if (longmessage) 
		{
		printf("\nPress any key to continue &");
		longmessage = FALSE;
		}
	else printf("\n&");
	gets(buf); 
	}