{
Module to Calculate Summary Statistics for Historical Fit based on:

Sterman, J.D., (1984) Appropriate Summary Statistics for Evaluating the 
Historical Fit of System Dynamics Models. Dynamica, 10 (Winter), 51-66

Rogelio Oliva.
September, 1995.
}

{Input Series}
input 1
	~	Units
	~	Rename this variable to the test data available in the .vdf file
		   The series has no value (no '=' sign) to indicate that it will be 
		   exogenous data.
	|

input 2
	~	Units
	~	Rename this variable to the test data available in the .vdf file
		   The series has no value (no '=' sign) to indicate that it will be 
		   exogenous data.
	|

historical := input 1
	~	Units
	~	Dummy variable to id series
	|

simulated := input 2
	~	Units
	~	Dummy variable to id series
	|

X  := simulated
	~	Units
	~		|

Y :RAW: := historical
	~	Units
	~		|

{Selection of points through the PICK function}
pick = IF THEN ELSE(Y = :NA:, 0, 1)
	~	Dimensionless
	~	Flag to id historical value available.
		   Takes a value of one for every data point available
	|

Xi = pick*X
	~	Units
	~	Simulated point entering calculations
	|

Yi = pick*Y
	~	Units
	~	Historical point entering calculations
	|

{Basic Accumulations}
count = INTEG(pick/dt,0)
	~	Dimensionless
	~	Counter for # of points
	|

Sum Xi = INTEG(Xi/dt,0)
	~	Units
	~	Sum of x's (simulated)
	|

SumX2 = INTEG(Xi*Xi/dt,0)
	~	Units*Units
	~	Sum of x^2 (simulated)
	|

Sum Yi = INTEG(Yi/dt,0)
	~	Units
	~	Sum of y's (historical)
	|

SumY2 = INTEG(Yi*Yi/dt,0)
	~	Units*Units
	~	Sum of y^2 (historical)
	|

SumXY = INTEG(Xi*Yi/dt,0)
	~	Units*Units
	~	Sum of x*y
	|

{Sum of Errors}
Sum XmY2 = INTEG((Xi-Yi)*(Xi-Yi)/dt,0)
	~	Units*Units
	~	Sum of Square Errors (x-y)^2
	|

Sum APE = INTEG(ABS(ZIDZ(Xi-Yi,Yi))/dt,0)
	~	Dimensionless
	~	Sum of Absolute Percent Errors
	|

Sum SPE = INTEG((ZIDZ(Xi-Yi,Yi)*ZIDZ(Xi-Yi,Yi))/dt,0)
	~	Dimensionless
	~	Sum of Square Percent Errors ((x-y)/y)^2
	|

{Calculation of Mean and Standard Deviation}
M X = ZIDZ(Sum Xi,count)
	~	Units
	~	Mean of x (sum x)/n
	|

M Y = ZIDZ(Sum Yi,count)
	~	Units
	~	Mean of y (sum y)/n
	|

MX2 = ZIDZ(SumX2,count)
	~	Units*Units
	~	Mean of x^2 (sum x^2)/n
	|

MY2 = ZIDZ(SumY2,count)
	~	Units*Units
	~	Mean of y^2 (sum y^2)/n
	|

Mxy = ZIDZ(SumXY,count)
	~	Units*Units
	~	Mean of x*y (sum x*y)/n
	|

Sx = SQRT(MX2-(M X*M X))
	~	Units
	~	Standard Deviation of x. Calculated using the 'hand computation' formula 
		   to calculate the standard deviation without prior knowledge of the mean. 
		   Sterman (1984), pg. 64
	|

Sy = SQRT(MY2-(M Y*M Y))
	~	Units
	~	Standard Deviation of y. Calculated using the 'hand computation' formula 
		   to calculate the standard deviation without prior knowledge of the mean. 
		   Sterman (1984), pg. 64
	|

r = ZIDZ(Mxy-(M X*M Y),Sx*Sy)
	~	Dimensionless
	~	Correlation coefficient. Calculated through the 'hand computation'. 
		   Sterman (1984) pg. 63
	|

{Decomposition of the Mean Square Error}
dif mea = (M X-M Y)*(M X-M Y)
	~	Units*Units
	~	Difference of Means (bias)
	|

dif var = (Sx-Sy)*(Sx-Sy)
	~	Units*Units
	~	Difference of variances
	|

dif cov = 2*Sx*Sy*(1-r)
	~	Units*Units
	~	Difference of covariances
	|

mse = dif mea + dif var + dif cov
	~	Units*Units
	~	Mean Square Error. The addition of the three components
	|

{Summary Statistics}
RMSPE = SQRT(ZIDZ(Sum SPE,count))
	~	Units
	~	Root Mean Square Percent Error
	~	:SUPPLEMENTARY 
	|

mape = ZIDZ(Sum APE,count)
	~	Dimensionless
	~	Mean Absolute Percent Error
	~	:SUPPLEMENTARY 
	|

rmse = SQRT(mse)
	~	Units
	~	Root Mean Square Error
	~	:SUPPLEMENTARY 
	|

um = ZIDZ(dif mea,mse)
	~	Dimensionless
	~	Bias inequality proportion
	~	:SUPPLEMENTARY 
	|

us = ZIDZ(dif var,mse)
	~	Dimensionless
	~	Variance inequality proportion
	~	:SUPPLEMENTARY 
	|

uc = ZIDZ(dif cov,mse)
	~	Dimensionless
	~	Covariance inequality proportion
	~	:SUPPLEMENTARY 
	|

r2 = r*r
	~	Dimensionless
	~	Correlation coefficient squared
	~	:SUPPLEMENTARY 
	|

dt = TIME STEP
	~	Time Units
	~		|

residuals = Xi-Yi
	~	Units
	~	Errors
	|

********************************************************
	.Control
********************************************************~
		Simulation Control Parameters
	|

FINAL TIME  = 55
	~	Time Units
	~	The final time for the simulation.
	|

INITIAL TIME  = 0
	~	Time Units
	~	The initial time for the simulation.
	|

SAVEPER  = 
        TIME STEP
	~	Time Units
	~	The frequency with which output is stored.
	|

TIME STEP  = 1
	~	Time Units
	~	The time step for the simulation.
	|

\\\---/// Sketch information - do not modify anything except names
V300  Do not put anything below this section - it will be ignored
*View 1
$0,0,Helvetica|14|B|0-0-0|0-0-0|0-0-0|-1--1--1|-1--1--1
12,1,0,429,257,130,23,8,4,0,0,0,0,0,0
Rename the inputs to the series you want to test and supply them in a .vdf file.
12,2,0,426,119,115,31,8,4,0,0,0,0,0,0
Ensure that time bounds cover the whole range of data availability; + 1 dt on the end
10,3,historical,351,191,42,12,0,3,0,0,-1,0,0,0
10,4,simulated,531,191,45,12,0,3,0,0,-1,0,0,0
10,5,input 1,171,192,40,20,8,3,0,2,0,0,0,0,-1--1--1,0-0-0,|0||255-128-128
10,6,input 2,673,189,40,20,8,3,0,2,0,0,0,0,-1--1--1,0-0-0,|0||255-128-128
1,7,5,3,0,0,0,0,0,0,0,-1--1--1,,1|(253,191)|
1,8,6,4,0,0,0,0,0,0,0,-1--1--1,,1|(611,189)|
///---\\\
:L<%^E!@
15:0,0,0,0
19:100,0
5:count