vqtam.m

MATLAB implementation of time series prediction Based on the VQTAM method described in the following

% MATLAB implementation of time series prediction
% Based on the VQTAM method described in the following papers:
% 
% G. A. Barreto & A. F. R. Araujo (2004)
% "Identification and Control of Dynamical Systems Using the Self-Organizing Map"
% IEEE Transactions on Neural Networks, vol. 15, no. 5.
%
%
% Authors: Guilherme A. Barreto & Leonardo Aguayo
% Date: August 30 2005

clear; clc; close all;

% Load data
load train_series_data.dat;    % Load the time series to be clustered

%% Building the input vectors from an univariate time series
p=5;       % Dimension of the input regression vector (length of the time window)
lap=p-1;   % Amount of overlapping between consecutive input vectors
Dw=buffer(train_series_data,p,lap);   % Build the data vectors
clear train_series_data
if lap>0,
     Dw=Dw(:,p:end)';  % Eliminate the first 'p-1' vectors with zeros)
else Dw=Dw';
end
Dw=fliplr(Dw);

Ynext=Dw(2:end,1);        % Target (one-step-ahead) values 
Dw=[Ynext Dw(1:end-1,:)]; % Input regressor vectors for Ynext

Dw=Dw+0.01*randn(size(Dw));  % Add some gaussian noise to the data

%---------------------------------------------------------
% SOM initialization and training
Mx = 16;            % Number of neurons in the X-dimension
My = 1;              % Number of neurons in the Y-dimension
msize = [Mx My];     % Size of 2-D SOM map
MASK=[0; ones(p,1)];   
sMap = som_randinit(Dw, 'msize', msize,'lattice','rect','shape','sheet');
sMap = som_seqtrain(sMap,Dw,'mask',MASK,'radius',[floor(0.5*max(msize)) 0],'sample_order','random','neigh','gaussian','trainlen',50);

%--------------- Organize the testing data ------------
load test_series_data.dat  % load target time series
Dw=buffer(test_series_data,p,lap); % Build input data vectors

if lap>0,
     Dw=Dw(:,p:end)';  % get rid of the first 'p-1' vectors with zeros
else Dw=Dw';
end
Dw=fliplr(Dw);

Ynext=Dw(2:end,1);        % Target (one-step-ahead) values 
Dw=[Ynext Dw(1:end-1,:)]; % Input regressor vectors for Ynext

Dw=Dw+0.01*randn(size(Dw));  % Add some gaussian noise to the data

% Find the winning neurons for testing set
Winners = som_bmus(sMap,Dw,1,MASK);

% Compute the predictions
Yhat = sMap.codebook(Winners,1); % We take only the Wout part of the weight vector

% Plot target and predicted time series
figure; plot(Ynext,'r-'); hold on; plot(Yhat,'b-'); hold off

% Normalized mean squared error
error=Ynext-Yhat;  % Prediction error
NMSE=var(error)/var(Ynext)