%% matlab cleaning
close all  
clear all
clc

%% data
% model
A = [0.2 0.3]; B = 0.5; N = 100;
theta0 = [-A B]'; % parameter vector
sigma_e = 1;
Ufin = 10;

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

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

%% simulation
% initialiazzation
Y = Yo; Yu = Yo; Yoe = Yo;
U = [Uo; U]; E = [Eo; E]; Yp= [NaN; NaN];
T = 1:(N+2);
% 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


%% 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(T,Y,'.', T,Yp,'.r'); axis([1 N+2 0 7]); 
        legend({'output','prediction'}),title('ARX model');
        xlabel('time [min]'), ylabel('Y_p(k)')
    subplot(3,2,2),    
        plot(T,Yoe,'.', T,Yu,'.r'); axis([1 N+2 0 7]); 
        legend({'output','prediction'}),title('OE model')
        xlabel('time [min]'), ylabel('Y_P(k) = Y(k|k-1)')                
    subplot(3,2,3),    
        plot(T,U,'.'); title('eXogeneous input')
        xlabel('time [min]'), ylabel('U(k)')
        axis([1 N+2 0 Ufin*1.1]); 
   subplot(3,2,4),    
        plot(T,E,'.'); title('Error signal')
        xlabel('time [min]'), ylabel('E(k)')
        axis([1 N+2 0 max(E)*1.1]);        
   subplot(3,1,3),    
        plot(T,Ep_oe,'O',T,Ep_arx,'x',T,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)')     