offline ( Gauss Newtone approach ) Prediction error method for ARMAX Model, Algorithm Implemented but not good results

% Simulation of an ARMAX model y(t)-0.7y(t-1)=0.5u(t-1)+e(t)+0.6e(t)

np=251;

a0=[1 -0.7 ]; b0=[0 0.5 ]; c0=[1 0.6]; % Set parameter values

u=randn(np,1)*sqrt(1); % Generate white noise with variance =1

e=randn(np,1)*sqrt(1); % Generate white noise with variance = 1

y=filter(b0,a0,u)+filter(c0,a0,e); % generate output data : y=(B/A)*u+(C/A)*e

z=[y u]; % store our inpute output data

x=rpem(z,[1 1 1 0 0 1],'ff',1); % build-in predction error method

x(end,:);

%% %%

theta=[-0.6 ;0.4 ; 0.6];% initiliazing

grad=[]; % dummy

iter_max=10;

for iter=1:iter_max

theta1=theta(1);

theta2=theta(2);

theta3=theta(3);

pe=filter([1 theta1],[1 theta3],y)-filter([0 theta2],[1 theta3],u);% Predictive errors, for the armax is : pe=(A/C)y(t)-(B/C)u(t)

for i=2:251 % Gradient for armax is given by grad(t)=(1/C)(-y(t-1),u(t-1),pe(t-1)) .

grad(1,:)=[0 0 0];% dummy value at time 1, to be deleted later

grad(i,:)=filter(1,[1 theta3],[-y(i-1) u(i-1) pe(i-1)]);% we started at i=2 so we can calculate using past value for other data

grad(1,:)=[];% delete first dummy row , we get 250 data

end

pe(1,:)=[];% delete first row , to get 250 data so we can compare it with grad

P=(grad'*grad)\(grad'*pe); % Gauss-Newtone Formula

theta=theta+0.8*P;%GaussNewtone Iterations

end

Answer 1

[Aman]

Hi Reda,

As per my understanding, you are using the Gauss-Newton method of optimizing the prediction error, but the result is not converging.

The matrix that has been created is near singular, so in that case the delta is not getting computed properly. In this case, you need to add a damping factor so that the proper delta is calculated. Refer to the below code, where I have added a damping factor of "1e-3".

%  Simulation of an ARMAX model  y(t)-0.7y(t-1)=0.5u(t-1)+e(t)+0.6e(t)
np=251;
a0=[1 -0.7 ];
b0=[0 0.5 ];
c0=[1 0.6]; % Set parameter values
u=randn(np,1)*sqrt(1); %  Generate white noise  with variance = 1
e=randn(np,1)*sqrt(1); % Generate white noise with variance = 1
y=filter(b0,a0,u)+filter(c0,a0,e); %  generate output data : y=(B/A)*u+(C/A)*e
z=[y u]; % store our input output data
%% 
y_prev = [0; y(1:end-1)]; % y(t-1)
u_prev = [0; u(1:end-1)]; % u(t-1)
e_prev = zeros(np, 1);
theta=[-0.6 ;0.4 ; 0.6];% initiliazing
iter_max=10;
error = [];
for iter=1:iter_max
    residual = y - (theta(1)*y_prev + theta(2)*u_prev + theta(3)*e_prev);
    J = [-y_prev, u_prev, e_prev];
    lambda = 1e-3; % Damping factor
    delta = ((J'*J + lambda*eye(size(J,2))) \ J') * residual;
    theta = theta + delta;
    e_prev = [0; residual(1:end-1)];
end

I hope it helps!