// P21EstReg1.sce
// Estimation 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]')

I_dataSim=2;                       // get data: 1=simulated, 2=real
ordE=2;                            // model order (valid for real data)

select I_dataSim
case 1 
  load _data/dataT11.dat Sim       // load reults of task T11
  yt=Sim.Cy.yt;                    // or use your own data yt, ut
  ut=Sim.Cy.ut;                    //  and define order of the model
  ord=Sim.Cy.ord;
case 2
  load _data/dataReg1.dat yt ut    // from file _data/dataReg1.dat
end

nd=length(yt);
if I_dataSim==1,    
  ordE=Sim.Cy.ord;
  th=Sim.Cy.th;
end
nPsi=2*ordE+3;                     // length of regression vector
  
V=1e-8*eye(nPsi,nPsi);             // comput. of information matrix
for t=(ordE+1):nd
  Ps=[yt(t) genpsi(t,ordE,yt,ut)]';
  V=V+Ps*Ps';
end
thE=inv(V(2:$,2:$))*V(2:$,1);      // point estimates of parameters
Est.Cy.ord=ordE;
Est.Cy.V=V;
Est.Cy.th=thE;

// results
if I_dataSim==1
  disp('Simulated parmeters')
  disp(th)
end
disp('Estimated parmeters')
disp(thE')

if I_dataSim==1
  if length(th)==length(thE)
    bar([th' thE])
    title('Simulated and estimated parameters')
  end
else
  bar(thE)
  title('Estimated parameters')
end

save _data/dataT21.dat Est
