%DMF Handout February 10, 1999

%yielddemo.m
%Demonstrate polynomial fit and simple regression using a multi-layer perceptron
%This demo is derived from the Netlab demo DEMMLP1

%
%	Neural Net: Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997)
%  Yield Curve:Copyright (c) Claudia Perlich, Andreas S Weigend (1999)

%  cd /home/sabai/aweigend/public_html/DMFS99/Notes/Class07
load crspdat.mat

clf
plot(xcrsp,tcrsp,'.');
title('December 1997 data from CRSP');
xlabel('Maturity in days (logarihm)')
ylabel('Percent yield')

disp('The figure diplays the true data (blue dots)')
disp('Press any key to continue.')
pause

hold on
plot(sort(xcrsp), yieldnet(sort(xcrsp)), '-r', 'LineWidth', 2) %plot the approx
hold off
disp('The figure diplays the true data (blue dots) and a neural net fit.')

%From now on we actually don't work with the original data but with 
%realization of the fit
%This allows you to build up intuitions about the effect of the number 
%of data points and the noise level

ndata = 10;		   	% Number of data points.
noise = 0.2;			% Standard deviation of noise distribution.

randn('state', 1);   %initialize the random number generator
x = [2.5:6/(ndata -1):8.5]';  %space x-values evenly
%x = 2.5 + 6 * rand(ndata,1)  %space x-values randomly

t=yieldnet(x) + noise*randn(length(x), 1); %generate target data with noise

%_____________________________Polynomial________________________________
disp(' ')
disp('Now we do polynomial fits, starting with a linear model.')
disp('Press any key to continue.')
pause
plotvals = [2.5:0.01:8.5]';

for orderpoly = 1:length(x)
   [p,s] = polyfit(x,t,orderpoly);
   ypoly = polyval(p,plotvals);
   plot(x,t,'.'); hold on; %plot the points
   plot(plotvals, ypoly, ':m', 'LineWidth', 2)
   title(orderpoly)
   xlabel('Maturity in days (logarihm)')
   ylabel('Percent yield (Polynomial fit)')
   pause; clf
end

%___________________________Neural Network_______________________________
% Set up network parameters.
disp('Setting up network...')
nin = 1;					% Number of inputs.
nhidden = 3;			% Number of hidden units.
nout = 1;				% Number of outputs.
alpha = 0.01;			% Coefficient of weight-decay prior, typical value 0.01

% Create and initialize network weight vector.

net = mlp(nin, nhidden, nout, 'linear', alpha);

% Set up vector of options for the optimiser.

options = zeros(1,18);
options(1) = 1;			% This provides display of error values.
options(14) = 120;		% Number of training epochs

disp(['The network has ', num2str(nhidden), ' hidden units and a weight decay'])
disp(['coefficient of ', num2str(alpha)])
disp('After initializing the network, we train it use the scaled conjugate')
disp(['gradients algorithm for ', num2str(options(14)), ' epochs.'])
disp(' ')
disp('Press any key to continue')
pause

% Train using scaled conjugate gradients (scg)
[net, options] = netopt(net, options, x, t, 'scg');

disp(' ')
disp('Now we plot the data, underlying function, and network outputs')
disp('on a single graph to compare the results.')
disp(' ')
pause

% Plot the data, the original function, and the trained network function.
y = mlpfwd(net, plotvals);
fh1 = figure;
plot(x, t, 'ob')
hold on
%axis([0 1 -1.5 1.5])
%[fx, fy] = fplot('sin(2*pi*x)', [0 1]);
plot(plotvals, yieldnet(plotvals), '-.r', 'LineWidth', 2)
plot(plotvals, y, '-k', 'LineWidth', 2)
plot(plotvals, ypoly, ':m', 'LineWidth', 2)

legend('data (realizations)', 'true function', ...
  'network trained to data','polynomial fit');
title('Yield curve example');
xlabel('Maturity in days (logarihm)')
ylabel('Percent yield')
hold off

disp('Demo done... returning the command window to the user...')