clear all
close all
clc

load heater

y1=ye-mean(ye);
u1=ue-mean(ue);
y2=yv-mean(yv);
u2=uv-mean(uv);

m=arx([y1' u1'],[2 2 4]);

[num,den]=tfdata(m,'v');
th=[-den(2:3) num(5:6)];

ys=zeros(size(y2));
for t=5:length(y2)
    % w=[y2(t);y2(t-1);u2(t-2);u2(t-3)]; % 1-step prediction
    w=[ys(t);ys(t-1);u2(t-3);u2(t-4)]; % simulation
    ys(t+1)=th*w;
end

% ys=sim(m,u2')';

RMSE=norm(y2(101:1000)-ys(101:1000))/sqrt(900)

figure
subplot(2,1,1),hold on,box on
plot(y2)
plot(ys,'r')
xlim([0 500])
subplot(2,1,2)
plot(u2)
xlim([0 500])

