%% matlab cleaning
close all  
clear all

%% model
A = [0.2 0.3]; B = 0.5; N = 100;
theta0 = [-A B]'; % parameter vector

%% input signal
E = randn(N,1); U = 10*ones(N,1);

%% initialization parameter
Yo = [0 0]'; Eo = [NaN NaN]'; Uo = [NaN 0]';

%% simulation
% initialiazzation
Y = Yo; Yu = Yo; Yoe = Yo;
U = [Uo; U]; E = [Eo; E]; Yp= [NaN; NaN];
% simulation cycle
for i = 2+(1:N)
    % modello arx
    phi = [Y(i-1) Y(i-2) U(i-1)]; % regressors are a selection of Z = subset of {Y,Y}
    Y(i,1) = phi*theta0 + E(i);      % ARX output
    Yp(i,1) = phi*theta0;            % ARX one-step prediction
    
    % modello oe
    phi_oe = [Yu(i-1) Yu(i-2) U(i-1)]; % regressors are not included in Z = subset of {Y,Y}
    Yu(i,1) = phi_oe*theta0;           % OE one-step prediction
    Yoe(i,1) = Yu(i)+E(i);             % OE output
end
% additional information removal
Y(1:2)=[];U(1:2)=[];E(1:2)=[]; Yp(1:2) = [];
Yu(1:2)=[]; Yoe(1:2)=[]; 
%% precition error evaluation
Ep_arx = Y - Yp;
Ep_oe = Yoe - Yu;
%% data save
save ('Es1Outcome' ,'Y','Yoe','U', 'E','theta0')

%% draw figures
figure (1)
    subplot(3,2,1), % define an area internal to the figure
        plot(1:N,Y,'.', 1:N,Yp,'.r'); axis([1 N 0 7]); 
        legend({'output','prediction'}),title('ARX model');
        xlabel('time [min]'), ylabel('Y_p(k)')
    subplot(3,2,2),    
        plot(1:N,Yoe,'.', 1:N,Yu,'.r'); axis([1 N 0 7]); 
        legend({'output','prediction'}),title('OE model')
        xlabel('time [min]'), ylabel('Y_P(k) = Y(k|k-1)')        
    subplot(3,2,3),    
        plot(1:N,U,'.'); title('eXogeneous input')
        xlabel('time [min]'), ylabel('U(k)')
   subplot(3,2,4),    
        plot(1:N,E,'.'); title('Error signal')
        xlabel('time [min]'), ylabel('E(k)')
        
   subplot(3,1,3),    
        plot(1:N,Ep_oe,'O',1:N,Ep_arx,'x',1:N,E,'.');         
        legend({'oe','arx','Simulation Error'});
        title('Prediction Error (E_p(k) = Y(k|k-1) -Y(k))')
        xlabel('time [min]'), ylabel('E_p(k)')     