%% Kittiphan Techakittiroj (27 Jan 1996)

%% Input
k = 0:1:100;
C = [50*sin(0.1*k)] + 700;
Sold = [ zeros(11,1); 100*ones(30,1); zeros(30,1); 50*ones(10,1); zeros(20,1)] + 1000;
CSold = C.*Sold';  CSold = CSold(:);        %% Make a new input from a non-linear operation
xa = [100*sin(0.1*k+3)] + 300;  xa = xa(:);
xb = [100*sin(0.2*k+3)] + 275;  xb = xb(:);
xc = [100*sin(0.3*k+3)] + 300;  xc = xc(:);
xd = [100*sin(0.4*k+3)] + 275;  xd = xd(:);
input = [ CSold xa xb xc xd];

%% Model
B0 = [   1    0    0    0    0];   A0 = [ 1];  %% delay 0
B1 = [   0    0    0    0    0];   A1 = [-1];  %% delay 1
B2 = [   0    0    0    0    0];               %% delay 2
B3 = [   0    0 -600    0    0];               %% delay 3
B4 = [   0    0    0    0 -550];               %% delay 4
B5 = [   0 -500    0    0    0];               %% delay 5
B6 = [   0    0    0    0    0];               %% delay 6
B7 = [   0    0    0 -400    0];               %% delay 7
AA = [ A0 A1];
BB = [ B0 B1 B2 B3 B4 B5 B6 B7];
Model = arx2th( AA, BB, 1, 5);
y = idsim(input, Model);

%% ARX
B0 = [ 0.3307 0    0    0    0];   A0 = [    1];  %% delay 0
B1 = [   0    0    0    0    0];   A1 = [-0.7850];  %% delay 1
B2 = [   0    0    0    0    0];               %% delay 2
B3 = [   0    0 3425    0    0];               %% delay 3
B4 = [   0    0    0    0 2558];               %% delay 4
B5 = [   0 1107    0    0    0];               %% delay 5
B6 = [   0    0    0    0    0];               %% delay 6
B7 = [   0    0    0 969.4   0];               %% delay 7
AA = [ A0 A1];
BB = [ B0 B1 B2 B3 B4 B5 B6 B7];
Model = arx2th( AA, BB, 1, 5);
y1 = idsim(input, Model);

%% RARX
B0 = [ 8.8454 0    0    0    0];   A0 = [    1];  %% delay 0
B1 = [   0    0    0    0    0];   A1 = [-0.62];  %% delay 1
B2 = [   0    0    0    0    0];               %% delay 2
B3 = [   0    0 7381    0    0];               %% delay 3
B4 = [   0    0    0    0 -5871];              %% delay 4
B5 = [   0 8197    0    0    0];               %% delay 5
B6 = [   0    0    0    0    0];               %% delay 6
B7 = [   0    0    0 -9119   0];               %% delay 7
AA = [ A0 A1];
BB = [ B0 B1 B2 B3 B4 B5 B6 B7];
Model = arx2th( AA, BB, 1, 5);
y2 = idsim(input, Model);

%% Plot Difference
plot(k, y, 'y', k, y1, 'r', k, y2, 'g');
title('an example output from the simulation');
ylabel('profit ($)');  xlabel('time (day)');  grid;
