%Retirar os valores

x = datasetATDPL1;

data_labels= x(2:length(x), 1); %datas
valores = [x{2:length(x), 2}]; %valores referentes a cada data

N = length(valores);
t = (0: N-1);
tt = (0:2*N-1); %escala temporal para a previsao

%figure(1)
%plot (t, valores, '-+')

haNaN = any(isnan(valores));  %Ha NaN por colunas?
colNaN = find(isnan(valores)); % Colunas com NaN

%Elimina linhas com NaN e reconstroi
%Sugestao:
% Reconstruir as linhas eliminadas usando extrapolacao

valores_1 = valores;

if any(isnan(valores))
    ind = find(isnan(valores));
    for k = 1 : length(ind)
        tt = t(ind(k) - 4 : ind(k) - 1);
        xx = valores_1(ind(k) - 4 : ind(k) - 1);
        aux = interp1(tt,xx,t(ind(k)),'pchip','extrap');
        if aux < 0
            valores_1(ind(k)) = interp1(tt,xx,t(ind(k)),'linear','extrap');
        else
            valores_1(ind(k)) = aux;
        end
    end
end

%figure(1)
%plot (t, valores_1, '-+')


%Media e desvio padrao de cada grupo de 31
mu = zeros(12,1);
aux = 1;
auxMax = 31;

for k = 1 : 12
    mu(k,1) = mean(valores_1(aux:auxMax));
    auxMax = min(auxMax + 31, length(valores_1));
    aux = aux + 31;

end

sigma = zeros(12,1);
aux_1 = 1;
auxMax_1 = 31;
for k = 1 : 12
    vetor_aux = valores_1(aux_1:auxMax_1);
    sigma(k,1) = std(vetor_aux);
    auxMax_1 = min(auxMax_1 + 31, length(valores_1));
    aux_1 = aux_1 + 31;

end

%outliers

aux_2 = 1;
auxMax_2 = 31;
valores_s_outliers = valores_1; %Substituicao dos outliers

for k = 1 : 12
    vetor_aux = valores_1(aux_2:auxMax_2);
    N = length(vetor_aux);
    MeanMat = repmat(mu(k),N,1); %Repete a media (mu) em N linhas
    SigmaMat = repmat(sigma(k),N,1); %Repete a media (sigma) em N linhas
    outliers = find(abs(vetor_aux - MeanMat) > 3 * SigmaMat); %identifica os outliers
    nout = length(outliers); %numero de outliers
    if nout
        for x = aux_2: auxMax_2
            disp(k);

            if valores_s_outliers(x) - mu(k) > 3* sigma(k)
                disp(k);
               valores_s_outliers(x) = mu(k) + 2.5*sigma(k);
            elseif mu(k) - valores_s_outliers(x) > 3*sigma(k)
               valores_s_outliers(x) = mu(k) - 2.5*sigma(k);
           end
        end
    end
    auxMax_2 = min(auxMax_2 + 31, length(valores_1));
    aux_2 = aux_2 + 31;
end

%disp(valores_s_outliers);
%figure(1)
%plot (t, valores_1, '-+',t,valores_s_outliers,'-*')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ficha 4 e pergunta 2 do mini projeto

%2.2 e 2.3 MINIPROJETO

%---Estimar tendencia parametrica (polinomial)
%trend de ordem 0
valores_s_outliers_t0 = detrend(valores_s_outliers,'constant'); %media da serie
tr1_0 = valores_s_outliers - valores_s_outliers_t0;
%trend de ordem 1
valores_s_outliers_t1 = detrend(valores_s_outliers,'linear'); %serie sem a tendencia
tr1_1 = valores_s_outliers - valores_s_outliers_t1;

%2.4 Rpresentar graficamente

%figure(2) %Para ver mudar numero
%subplot(211)
%plot(t,valores_s_outliers,'-+',t,tr1_1,'-',t,valores_s_outliers_t1,'-');
%legend('Serie 1 regularizada', 'componente da tendência','Série temporal sem a tendência de grau 1');
t = t';
%vai se usar o t' e nao o t para que t e valores_s_outliers tenham as mesmas dimensoes
p = polyfit(t',valores_s_outliers,1); %coeficiente dos polinomios
tr1_2 = polyval(p,t); %valores resultantes
valores_s_tend = valores_s_outliers - tr1_1;

%2.5 Representacao grafica

%figure(1)
%subplot(211)
%plot(t,valores_s_outliers,'-+',t,tr1_1,'-*');
%title('Serie 1 (+) e trend(*) de grau 2')
%xlabel('t [h]');

%subplot(212)
%plot(t,valores_s_tend,'-o');
%title('Serie 1 sem (o) trend de grau 2')
%xlabel('t [h]');

%2.6 e %2.7
%N = length(valores);
%Estimar sazonalidade

h0 = repmat((1:91)', 4, 1); %trimestre
%h0 = cat(1, h0, h0(1:4));
sX = dummyvar(h0); %gera matriz dummy, matriz de 0's e 1's
Bs1 = sX\valores_s_outliers_t1'; %Obtem valores para cada trimestre %valores_s_outliers_t1??????
st1 = sX * Bs1; %componente sazonal

valores_s_sazonalidade = valores_s_outliers' - st1;

%2.8 Representacao grafica

%figure(10)
%subplot(211)
%plot(t,valores_s_outliers,'-+',t,valores_s_sazonalidade','-o');
%subplot(212)
%plot(t,st1,'-+');
%title('Sazonalidade da Serie 1')
%xlabel('t [h]');

%2.9

%-- Estimar a componente irregular
irreg = valores_s_outliers' - tr1_2 - st1;
%--- serie temporal sem a componente irregular
valores_s_irreg = valores_s_outliers' - irreg;

%2.10 Representacao grafica

%figure(1)
%subplot(211)
%plot(t,valores_s_outliers,'-+',t,valores_s_irreg','-o');
%title('Serie 1 com (+) e sem (o) componente irregular')
%subplot(212)
%plot(t,irreg,'-*');
%title('Componente irregular da serie 1')
%xlabel('t [h]');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ficha 5 e pergunta 3 do mini projeto

%3.1
%3.2

adftest(valores_s_outliers) %teste de estacionaridade da serie regularizada
adftest(st1) %teste de estacionaridade da componente sazonal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.3%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure(1); %representacao do FAC e do FACP

y1 = st1(1:91);

%subplot(211)
%autocorr(y1) %indicador para a ordem de nc (modelo MA)

%subplot(212)
%parcorr(y1) %indicador para a ordem de na (modelo RA)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.4%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

id_y1 = iddata(y1,[],1,'TimeUnit','day');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.5%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opt1_AR = arOptions('Approach','ls'); %Opçoes do modelo AR
na1_AR = 30; %historico da variavel
model1_AR = ar(id_y1,na1_AR,opt1_AR); %modelo AR
pcoef1_AR = polydata(model1_AR); %parametros do modelo AR

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.6%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Simulaçao do modelo AR
y1_AR = y1(1:na1_AR);
for k = na1_AR + 1 : 91
    y1_AR(k) = sum(-pcoef1_AR(2:end)'.*flip(y1_AR(k-na1_AR : k-1)));
end
y1_AR2 = repmat(y1_AR,4,1);

%Simulacao do modelo AR com forecast
y1_ARf = forecast(model1_AR,y1(1:na1_AR),91-na1_AR);
y1_ARf2 = repmat([y1(1:na1_AR); y1_ARf],4,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.7%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Representacoes graficas

%t=(1:182)';
%st1_=st1(1:182);

%figure(1) %compara a componente sazonal com a sua estimacao
%plot(t,st1,'+-',1:size(y1_AR2+tr1_2),y1_AR2+tr1_2,'-o');

%figure(2) %compara a componente sazonal com a sua estimacao
%plot(t,st1,'+-',t,y1_AR2,'-o',t,y1_ARf2,'-*');
%xlabel('t [h]');
%title('Componente sazonal 1 (-+) e Estimacao com o model AR (-o)');

%figure(3)
%plot(t,valores_s_outliers,'-+',t,y1_AR2+tr1_2,'-o');
%xlabel('t [h]');
%title('Serie 1 (-+) e estimaçao com o modelo AR (-o)');

%metrica para analise
E1_AR = sum((valores_s_outliers-(y1_AR2+tr1_2)).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.8%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tt = (0:364*2-1);

%tentativa de encontrar uma tendencia cos()
t1_prev = polyval(p, tt); %tendencia adaptada para 730 (linear) %calcula a endencia para 2N
media_G = mean(valores_s_outliers);
valores_normal = valores_s_outliers - media_G;
%plot(t, valores_s_outliers, '-+', t, valores_normal, 'r-+');
amplitude = 20;
periodo = 235;
acerto = -20;
tend_cos = amplitude * cos((2*pi/periodo)*tt - acerto) + t1_prev;
%plot(tt, tend_cos, '-+', t, valores_s_outliers, 'r-+');

figure(4) %Faz a previsao para 2N
aux3 = repmat(y1_AR2,2,1);
plot(t,valores_s_outliers,'-+',tt,aux3+tend_cos', '-o');
xlabel('t [h]');
title('Serie 1 (-+) e Previsao com o modelo AR (-o)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.9%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

opt1_ARMAX = armaxOptions('SearchMethod','auto'); %Opçoes do modelo ARMAX
na1_ARMA = 5; %numero de parcelas %mais valores
nc1_ARMA = 1; %ruido
model1_ARMA = armax(id_y1,[na1_ARMA nc1_ARMA],opt1_ARMAX);

[pa1_ARMA,pb1_ARMA,pc1_ARMA] = polydata(model1_ARMA);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.10%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

e = randn(91,1); %ruido branco

y1_ARMA = y1(1:na1_ARMA); %Simulacao do Modelo ARMA
for k = na1_ARMA + 1 : 91
    y1_ARMA(k) = sum(-pa1_ARMA(2:end)'.*flip(y1_ARMA(k - na1_ARMA : k - 1))) + sum(pc1_ARMA'.*flip(e(k - nc1_ARMA : k)));
end
y1_ARMA2 = repmat(y1_ARMA,4,1); %que valor aqui?????????

%Simulacao do Modelo ARMA com forecast

y1_ARMAf = forecast(model1_ARMA,y1(1:na1_ARMA), 91 - na1_ARMA); %que valor aqui?????????
y1_ARMAf2 = repmat([y1(1:na1_ARMA); y1_ARMAf],4,1); %que valor aqui?????????


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.11%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure(1) %compara a componente sazonal com a sua estimacao
%plot(t,st1,'+-',t,y1_ARMA2,'-o',t,y1_ARMAf2,'-*');
%xlabel('t [h]');
%title('Componente sazonal 1 (-+) e Estimacao com o model AR (-o)');

%figure(1)
%plot(t,valores_s_outliers,'-+',t,y1_ARMA2+tr1_2,'-o');
%xlabel('t [h]');
%title('Serie 1 (-+) e estimaçao com o modelo AR (-o)');

%metrica para analise
E1_AR = sum((valores_s_outliers-(y1_ARMA2+tr1_2)).^2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.12%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%provavelmente esta mal, ver o que e !!!!!!!!!!!!!

%figure(1) %Faz a previsao para 2N
%plot(t,valores_s_outliers,'-+',tt,repmat(y1_ARMA2,4,1) + tr1_2_2, '-o');
%xlabel('t [h]');
%title('Serie 1 (-+) e Previsao com o modelo ARMA (-o)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.13%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Estimcao de um modelo ARIMA

%Parametros do modelo ARIMA
p1_ARIMA = 3;  %?????
D1_ARIMA = 1; %?????
q1_ARIMA = 1; %?????

%Obtem a estrutura do modelo ARIMA
Md1 = arima(p1_ARIMA,D1_ARIMA,q1_ARIMA);

%Estimacao do modelo ARIMA ?????????????????
%EstMd1 = estimate(Md1,valores_s_outliers(1:N) , 'Y0', valores_s_outliers(1:p1_ARIMA +1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.14%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Simulacao do modelo ARIMA
%y1_ARIMA = simulate(EstMd1,N);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%3.15%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figure(1)
%plot(t,valores_s_outliers,'-+',t,y1_ARIMA,'-o');
%xlabel('t [h]');
%title('Serie 1 (-+) e Estimacao com o modelo ARIMA (-o)')

%Metrica para analise ARIMA
%E1_ARIMA = sum((valores_s_outliers - y1_ARIMA(1:N)).^2);