// P11SimReg1.sce
// Simulation of scalar regression model of order n    
// ------------------------------------------------
[u,t,n]=file();                         // find working directory
chdir(dirname(n(2)));                   // set working directory
clear("u","t","n")                      // clear auxiliary data
exec("ScIntro.sce",-1),mode(0)          // intro to sesion
deff('ps=genpsi(t,n,y,u)','ps=[y(t-(1:n)),u(t-(0:n)) 1]')

nd=100;                                 // number of data
I_typU=1;                               // type of input

// model
Sim.Cy.ord=2;                           // model order
Sim.Cy.tha=[.6 -.2 .1 -.5 .2];          // parameters at y
Sim.Cy.thb=[1 .5 -.3 -.1 .1 -.3];       // parameters at u
Sim.Cy.thk=-1;                          // model constant
Sim.Cy.cv=.01;                           // variance of noise

// model of input signal
select I_typU                           // choice of input
  case 1, ut=rand(1,nd,'u');            // random, uniform
  case 2, ut=ones(1,nd);                 // ones
  case 3                            // two steps
    ut=[ones(1,fix(nd/2)), -ones(1,fix(nd/2)+1)]; 
  case 4                                // random jumps
    ut(1)=1; j=1;
    for i=2:nd
      if rand(1,1,'u')>.85, j=rand(1,1,'n'); end
      ut=[ut j];
    end
  case 5, ut=abs(sin(10*(1:nd)/nd));    // sinus
end 
yt=zeros(1,nd);
ord=Sim.Cy.ord;
cv=Sim.Cy.cv;
th=[Sim.Cy.tha(1:ord) Sim.Cy.thb(1:(ord+1)) Sim.Cy.thk]; 
                                        // regression coefficients
Sim.Cy.th=th;

// time loop
for t=(ord+1):nd                        // time loop
  ps=genpsi(t,ord,yt,ut);               // regression vector
  yt(t)=ps*th'+cv*rand(1,1,'n');        // model simulation
end
Sim.Cy.yt=yt;
Sim.Cy.ut=ut;

// results
s=1:nd;                                 
plot(s,yt(s),s,ut(s))
legend('output','input');
title('Simulation with regression model')
set(gca(),'data_bounds',[1 nd min([yt,ut])-.1 max([yt,ut])+.1])

save _data/dataT11.dat Sim

