% testLinearAdditiveSDE: script to test the correctness of linearAdditive
%
% This script file runs the function linearAdditive to compute
%   the mean and variance of motion according to a linear ODE with
%   additive white noise:
%
%	dx = a * x * dt + b dW
%
% Constants a,b may be modified within the script.  It is also possible
%   to run a 2D version of the problem, which is just the 1D version
%   with variable x_1 + x_2 (the orthogonal direction x_1 - x_2 has
%   zero dynamics).
%
% The analytic mean and variance of such a flow is given by
%
%       E[x]   = exp(a * t) E[x(0)]
%	Var[x] = 0.5 / a * ((b^2 + 2 * a * Var[x(0)]) * exp(2 * a * t) - b^2)
%
%   which is taken from the end of a chapter of lecture notes at
%
% 	http://www.eng.tau.ac.il/~liptser/lectures/lect_new10.pdf
%
%   The script plots these analytic solutions for comparison purposes.
%
% Note that many problem parameters can be modified in the routine
%   linearAdditive.  However, modifications there may cause the
%   analytic solutions plotted here to be incorrect.

% Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell, 8/26/04

%---------------------------------------------------------------------------
% Problem parameters
a = 1.0;
b = 1.0;
tf = 0.25;
dim = 1;
accuracy = 'medium';

% Where to put the legend (according to help entry for legend)?
legendLocation = 4;

%---------------------------------------------------------------------------
% Compute the mean state at final time.
%   This call will also give us a grid structure to work with.
payoff = 'x';
[ dataX, g ] = linearAdditiveSDE(payoff, a, b, tf, dim, accuracy);

% We don't need the figure generated by linearAdditive.
close;

%---------------------------------------------------------------------------
% What is our dependent variable, 
%   and what grid nodes should we examine for results?
switch(dim)
 case 1
  % Use the only dimension we are given.
  indexX = 1 : g.N(1);
  stateX = g.xs{1};
 case 2
  % The operative dimension is the diagonal x_1 = x_2.
  indexX = 1 : g.N(1) + 1 : prod(g.N);
  stateX = (g.xs{1}(indexX) + g.xs{2}(indexX)) / sqrt(2);
 otherwise
  error('Dimension %d not currently allowed.', dim);
end

%---------------------------------------------------------------------------
% Analytic results.
mean_t0 = stateX;
var_t0 = zeros(size(mean_t0));

mean_tf = exp(a * tf) * mean_t0;
var_tf = 0.5 / a * ((b^2 + 2 * a * var_t0) * exp(2 * a * tf) - b^2);

%---------------------------------------------------------------------------
% Plot the results for the mean.
figure;
subplot(2,1,1);
h = plot(stateX, dataX(indexX), 'b-o', stateX, mean_tf, 'r--x');
legend('from PDE', 'analytic', legendLocation);
xlabel('x0');  ylabel('mean');
title([ 'Expected final state at t = ' num2str(tf) ]);

subplot(2,1,2);
h = [ h; plot(stateX, mean_tf - dataX(indexX), 'k-*') ];
xlabel('x0');  ylabel('error (analytic - calculated)');

%---------------------------------------------------------------------------
% Compute the variance at the final time.
payoff = 'x^2';
[ dataX2, g ] = linearAdditiveSDE(payoff, a, b, tf, dim, accuracy);

% We don't need the figure generated by linearAdditive.
close;

dataVar = dataX2 - dataX.^2;

%---------------------------------------------------------------------------
% Plot the results for the variance.
figure;
subplot(2,1,1);
h = [ h; plot(stateX, dataVar(indexX), 'b-o', stateX, var_tf, 'r--x') ];
legend('from PDE', 'analytic', legendLocation);
xlabel('x0');  ylabel('variance');
title([ 'Variance in state at t = ' num2str(tf) ]);

subplot(2,1,2);
h = [ h; plot(stateX, var_tf - dataVar(indexX), 'k-*') ];
xlabel('x0');  ylabel('error (analytic - calculated)');

% Figures for the documentation need thick lines and no symbols
%set(h, 'LineWidth', 2, 'Marker', 'none');

