

% nostri dati

load Dati_Matrice_norm.mat
Ts=8*60;

tic

for k=1:1%length(Dati) 
    
    y = Dati(k).ytemp_norm; % necessita dati a media 0. in part il fit.
    U = Dati(k).Udati_per_t_norm;   

           
    i_start= 50;  % 800 per ident, 600 per validazione. si scartano gli ultimi.
    i_stop = 850;
    Up=U(i_start:i_stop,:);     % p per passato
    yp=y(i_start:i_stop,:);

    Uf=U(i_stop+1:1450,:);       % f per futuro (per validazione) % qua tagliare la parte rumorosa finale.TO DO.
    yf=y(i_stop+1:1450,:);

    dat_id = iddata(yp,Up,Ts);
    dat_val = iddata(yf,Uf,Ts);



 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% STIMA PEM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% dat=iddata(y,U);
% dat_test=iddata(ytest_tot,Utest);


dat=dat_id;
dat_test=dat_val;
y = yp;

t=0;

n=800; % dati per id
p=30;
m=32;
mPEM=min(15,ceil((9/10)*(n+p)/(m+2)))%20;%max ordine PEM


for i1=1:mPEM;
    
    display(['iterazione' num2str(i1)])
            i2=i1*ones(1,m);
            i3=i1;
            t=t+1;
            Mpem=pem(dat,[i1 i2 i3 zeros(1,m+1) ones(1,m)],'LimitError',0);
            predd=predict(Mpem,dat);predd=iddata(predd);predd=predd.y;
            RM=sum((predd-y).^2);npar=max(size(Mpem.ParameterVector));ndat=Mpem.EstimationInfo.DataLength;
            L1(t)=ndat*log(RM/ndat)+log(ndat)*npar;%BIC
            L2(t)=ndat*log(RM/ndat)+[ndat*(ndat+npar)]/[ndat-npar-2];%AICC
            ll1(t,:)=[i1 i2 i3];

%             MpemOR=pem(dat,[i1 i2.*Ind(jj,(2:end)) i3 zeros(1,m+1) ones(1,m)],'LimitError',0);
%             predd=predict(MpemOR,dat);predd=iddata(predd);predd=predd.y;
%             RM=sum((predd-y).^2);npar=max(size(MpemOR.ParameterVector));ndat=MpemOR.EstimationInfo.DataLength;
%             L1b(t)=ndat*log(RM/ndat)+log(ndat)*npar;%BIC+1/2Or


    end
    [q,w1]=min(L1);ll1(w1,:)
    Mpem_BIC=pem(dat,[ll1(w1,:) zeros(1,m+1) ones(1,m)],'LimitError',0);

%     [q,w1]=min(L1b);i1=ll1(w1,1);i2=i1;i3=i1;ll1(w1,:)
%     Mpem=pem(dat,[i1 i2.*Ind(jj,(2:end)) i3 zeros(1,m+1) ones(1,m)],'LimitError',0);

    [q,w2]=min(L2);ll1(w2,:)
    Mpem_AICC=pem(dat,[ll1(w2,:) zeros(1,m+1) ones(1,m)],'LimitError',0);

    i1=ll1(w3,1);i2=i1*ones(1,m);i3=i1

%     MpemOR=pem(dat,[ll1(w3,1)  ll1(w3,2).*Ind(jj,(2:end))  ll1(w3,3)  zeros(1,m+1) ones(1,m)],'LimitError',0);



end

toc


%%


ystim = sim(Mpem,Uf);

compare(dat_val,Mpem,50);

figure
hold on
plot(yf,'r')
plot(ystim,'b')
hold off

