Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % 1)
- % 1.1) Leitura e representação dos resultados lidos
- data = importdata('dataset_ATD_PL7.csv');
- leituras = data.data;
- datas = data.textdata(2:366);
- N = length(leituras);
- t = 0:N-1;
- figure(1);
- plot(t, leituras, '-+');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- legend('Leituras ao longo do ano','Location', 'NorthEast');
- % 1.2) Verificação da existência de NaN através de extrapolação
- x1 = leituras;
- if any(isnan(x1))
- ind = find(isnan(x1));
- x1r = x1;
- for k=1: length(ind)
- tt = t(ind(k)-4 : ind(k)-1);
- xx = x1r (ind(k)-4 : ind(k)-1);
- x1r(ind(k)) = interp1(tt, xx, t(ind(k)), 'pchip', 'extrap');
- end
- end
- figure(2);
- plot(t, x1, '-+', t, x1r, '-o');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- legend('com NaN', 'sem NaN', 'Location', 'NorthEast');
- % 1.3) e 1.4) Verificação da existência de outliers através de uma janela deslizante
- % valores da média e desvio padrão são calculados a cada novo mês
- % substituição dos outliers por valores adequados e sua respetiva
- % representação
- % variáveis auxliares para a janela deslizante
- tam = N;
- passo=30;
- conta=1;
- mes=1;
- % inicializar vetor a 0
- s = size(x1r);
- x1ro = zeros(s);
- for i = 1: passo: N
- % caso em que o mês tem 30 dias, até ao mês 12
- if (mes <= 12)
- new_x1r = x1r(i:i+(passo-1));
- x1ro(i:i+(passo-1)) = new_x1r;
- % cálculo das médias do mês
- mu1 = mean(new_x1r);
- sigma1 = std(new_x1r);
- % verifica 1 a 1 cada dia do mes
- for j = 1: 1: 30
- if (abs(new_x1r(j) - mu1) > 3*sigma1)
- disp('outlier')
- if x1ro(i+j-1) > mu1
- disp(x1ro(i+j-1))
- x1ro(i+j-1) = mu1 + 2.5*sigma1;
- disp(x1ro(i+j-1))
- else
- disp(x1ro(i+j-1))
- x1ro(i+j-1) = mu1 - 2.5*sigma1;
- disp(x1ro(i+j-1))
- end
- end
- end
- mes = mes + 1;
- % caso do resto dos 5 dias que sobram para o total dos 365 dias
- else
- new_x1r = x1r(i:i+(5-1));
- x1ro(i:i+(5-1)) = new_x1r;
- % cálculo das médias do mês
- mu1 = mean(new_x1r);
- sigma1 = std(new_x1r);
- % verifica 1 a 1 cada dia do mes
- for j = 1: 1: 5
- if (abs(new_x1r(j) - mu1) > 3*sigma1)
- disp('outlier')
- if x1ro(i+j-1) > mu1
- disp(x1ro(i+j-1))
- x1ro(i+j-1) = mu1 + 2.5*sigma1;
- disp(x1ro(i+j-1))
- else
- disp(x1ro(i+j-1))
- x1ro(i+j-1) = mu1 - 2.5*sigma1;
- disp(x1ro(i+j-1))
- end
- end
- end
- mes = mes + 1;
- end
- end
- figure(3);
- plot(t, x1r, '-+', t, x1ro, '-o');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- legend('com outliers', 'outliers corrigidos', 'Location', 'NorthEast');
- %% 2
- N=length(x1ro);
- t=0:N-1;
- % 2.1) estimar serie temporal sem a componente de tendencia
- % aproximações polinomiais de grau 0 e 1, através da função detrend
- % 2.2) obter a componente
- x1ro_t0 = detrend(x1ro, 'constant');
- x1ro_t1 = detrend(x1ro, 'linear');
- % 2.3) Representar graficamente a série temporal regularizada, com
- % componente e sem componente de tendência
- % série temporal com tendência de grau 0
- figure(4)
- subplot(211)
- plot(t, x1ro, '-+', t, x1ro - x1ro_t0, '-*');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal com componente de tendência de grau 0');
- legend('Série temporal', 'Componente de tendência', 'Location', 'NorthEast');
- subplot(212)
- plot(t, x1ro_t0, '-o');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal sem componente de tendência de grau 0');
- legend('Sem componente de tendência', 'Location', 'NorthEast');
- % série temporal com tendência de grau 1
- figure(5)
- subplot(211)
- plot(t, x1ro, '-+', t, x1ro - x1ro_t1, '-*');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal com componente de tendência de grau 1');
- legend('Série temporal', 'Componente de tendência', 'Location', 'NorthEast');
- subplot(212)
- plot(t, x1ro_t1, '-o');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal sem componente de tendência de grau 1');
- legend('Sem componente de tendência', 'Location', 'NorthEast');
- % série com tendência de grau 5, caso polinomial
- % 2.4) obter a componente
- p = polyfit(t', x1ro, 5);
- tr1 = polyval(p, t');
- % 2.5) Representar graficamente a série temporal com e sem componente de
- % tendência
- figure(6)
- subplot(211)
- plot(t, x1ro, '-+', t, tr1, '-*');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal com componente de tendência de grau 5');
- legend('Série temporal', 'Com componente de tendência', 'Location', 'NorthEast');
- subplot(212)
- plot(t, x1ro-tr1, '-o');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal sem componente de tendência de grau 5');
- legend('Sem componente de tendência', 'Location', 'NorthEast');
- %{
- % SERIE SINUSOIDAL
- t_f = t.*24*3600;
- [Cm, tetam] = SerieFourier(t_f, x1ro, 360*24*60*60, 3);
- tr1 = 19.8784445884705*sin(2.23762861251355*10^(-17)*t);
- mt =[0 1];
- for k=1: length(mt)
- x1 = zeros(size(x));
- for m=0:mt(k)
- x1 = x1+Cm(m+1)*cos(m*2*pi/T0*t+tetam(m+1));
- end
- plot(t, x1, '-.b')
- end
- %}
- T0 = 365*24*60*60;
- m_max = 50;
- t_f = t.*24*3600;
- [Cm,tetam] = SerieFourier(t_f',x1ro,T0,m_max);
- mt =0:50;
- for k=1:length(mt)
- x1 = zeros(size(t));
- for m=0:mt(k)
- x1 = x1+ Cm(m+1)*cos(m*2*pi/T0*t_f+tetam(m+1));
- end
- end
- figure(40)
- plot(t_f,x1);
- hold on;
- tr1 = x1;
- %tr1 = 10*sin(t/1)+20;
- tr1 = tr1';
- figure(20)
- subplot(211)
- plot(t,tr1,'-',t,x1ro,'-*')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal com componente de tendência sinosuidal');
- legend('Série temporal', 'Componente de tendência', 'Location', 'NorthEast');
- subplot(212)
- plot(t, x1ro-tr1, '-o');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal sem componente de tendência sinosuidal');
- legend('Sem componente de tendência', 'Location', 'NorthEast');
- % 2.6) % Estimar a componente da sazonalidade MENSAL (30 dias)
- x1ro_t2 = x1ro-tr1;
- ho = (1:30)';
- s = size(ho)*13;
- s(2) = 1;
- x1ro_s = zeros(s);
- new_x1ro_t2 = zeros(s);
- new_x1ro_t2(1:365) = x1ro_t2(1:365);
- new_x1ro_t2(366:end) = x1ro_t2(1:25);
- % cálculo da componente da sazonalidade
- dv = repmat(ho, 13, 1);
- sx = dummyvar(dv);
- xtemp = sx\new_x1ro_t2;
- x1ro_s = sx * xtemp;
- x1ro_s = x1ro_s(1:365);
- % 2.7) Obter série sem a componente da sazonalidade
- % 2.8) Representar a série com e sem a componente da sazonalidade
- figure(7)
- subplot(211)
- plot(t, x1ro, '-+', t, x1ro-x1ro_s, '-*', t, tr1, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal com componente da sazonalidade');
- legend('Série temporal', 'Componente da sazonalidade', 'Location', 'NorthEast');
- subplot(212)
- plot(t, x1ro_s, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal sem componente da sazonalidade');
- legend('Sem da sazonalidade', 'Location', 'NorthEast');
- % 2.9) obter série temporal com e sem componente irregular
- irr_x1ro = x1ro - tr1 - x1ro_s;
- % 2.10) Representação da serie temporal com e sem componente irregular
- figure(8)
- subplot(211)
- plot(t, x1ro, '-+', t, x1ro-irr_x1ro, '-*')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal com componente irregular');
- legend('Com componente irregular', 'Location', 'NorthEast');
- subplot(212)
- plot(t, irr_x1ro, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Série temporal sem componente irregular');
- legend('Sem componente irregular', 'Location', 'NorthEast');
- %% 3)
- % 3.1) Obter componentes da série temporal resultantes da decomposição da
- % série temporal regularizada
- % 3.2) Verificar estacionaridade da série temporal através do adftest
- N = length(x1ro);
- t = (0:N-1)';
- tt = (0:2*N-1)';
- adftest(x1ro) % nao é estacionaria = 0
- adftest(x1ro_s) % sazonal deu 1- provavelmente é estacionaria
- % sazonal, as carateristicas repetem-se igualmente de 24 em 24 horas
- % componente sazonal da série
- y1 = x1ro_s(1:30);
- % 3.3) Representar graficamente FAC e FACP, através do autocorr e parcorr
- figure(9);
- subplot(211)
- autocorr(y1)
- subplot(212)
- parcorr(y1)
- % 3.4) Criar objeto iddata, com Ts=1
- id_y1 = iddata(y1, [], 1, 'TimeUnit', 'days');
- % 3.5) Modelo AR para a componente sazonal, resultado da FACP
- % definir na, usando arOptions, ar e polydata
- opt1_AR = arOptions('Approach', 'ls');
- na1_AR = 10; % repetir o processo ate encontrar numero aconselhavel, mais aceitavel
- model1_AR = ar(id_y1, na1_AR, opt1_AR)
- pcoef1_AR = polydata(model1_AR);
- y1_AR = y1(1:na1_AR);
- for k=na1_AR+1:30 % do ...(na1_AR+1) ao 30
- y1_AR(k) = sum(-pcoef1_AR(2:end)'.*flip(y1_AR(k-na1_AR: k-1)));
- end
- % 3.6) Teste gráfico usando parametros do modelo da componente sazonal e forecast
- y1_AR2 = repmat(y1_AR, 13, 1);
- y1_ARf = forecast(model1_AR, y1(1: na1_AR), 30-na1_AR);
- y1_ARf = repmat([y1(1:na1_AR); y1_ARf], 13, 1);
- % encurtar matrizes para 365 dias
- y1_AR2 = y1_AR2(1:365);
- y1_ARf = y1_ARf(1:365);
- figure(10)
- plot(t, x1ro_s, '-+', t, y1_AR2, '-o', t, y1_ARf, '-*');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Modelo AR');
- legend('Série temporal', 'AR calculado com for', 'AR calculado com forecast', 'Location', 'NorthEast');
- % 3.7) Validar o modelo comparando a série regularizada com o resultado da
- % combinação aditiva da estimação da componente sazonal + componente da
- % tendencia paramétrica da série. NOTA: mértica = erro^2 dos valores
- % medidos + erro^2 dos valores estimados
- % calcular erros (métrica)
- E1_AR = sum((x1ro_s-y1_AR2).^2);
- % 3.8) Gráfico do modelo AR para 2*tempo
- % aumentar para 2*365 dias
- %tr1_AR = polyval(p, tt);
- x1 = repmat(x1',2,1);
- tr1_AR = x1;
- %tr1_AR = 10*sin(tt/1)+20;
- st1_AR = repmat(y1_AR2, 2, 1);
- figure(11)
- plot(t, x1ro, '-+', tt, tr1_AR+st1_AR, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Modelo AR');
- legend('Série temporal', 'Previsão da série para o dobro da duração temporal', 'Location', 'NorthEast');
- % 3.9) Modelo ARMA para a componente sazonal, resultado da FAC
- % definir na, nc, usando armaxOptions, armax e polydata
- opt1_ARMAX = armaxOptions('SearchMethod', 'auto')
- ma1_ARMA = 5;
- mc1_ARMA = 1;
- model1_ARMA = armax(id_y1, [ma1_ARMA mc1_ARMA], opt1_ARMAX);
- [pa1_ARMA, pb1_ARMA, pc1_ARMA] = polydata(model1_ARMA);
- % ruido branco
- e = randn(30, 1);
- y1_ARMA = y1(1:ma1_ARMA);
- for k = ma1_ARMA + 1:30
- y1_ARMA(k) = sum(-pa1_ARMA(2:end)'.*flip(y1_ARMA(k - ma1_ARMA:k-1))) + sum(pc1_ARMA'.*flip(e(k-mc1_ARMA:k)));
- end
- % 3.10) Teste gráfico usando parametros do modelo da componente sazonal e forecast
- % feito o anterior automaticamente com a funcao forcast
- y1_ARMA2 = repmat(y1_ARMA, 13, 1);
- y1_ARMAf = forecast(model1_ARMA, y1(1:ma1_ARMA), 30-ma1_ARMA);
- y1_ARMAf2 = repmat([y1_ARMA(1:ma1_ARMA); y1_ARMAf], 13, 1);
- % encurtar matrizes para 365 dias
- y1_ARMA2 = y1_ARMA2(1:365);
- y1_ARMAf2 = y1_ARMAf2(1:365);
- figure(12);
- plot(t, x1ro_s, '-+', t, y1_ARMA2, '-o', t, y1_ARMAf2, '-*');
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Modelo ARMA');
- legend('Série temporal', 'ARMA calculado com for', 'ARMA calculado com forecast', 'Location', 'NorthEast');
- % 3.11) Validar o modelo comparando a série regularizada com o resultado da
- % combinação aditiva da estimação da componente sazonal + componente da
- % tendencia paramétrica da série. NOTA: mértica = erro^2 dos valores
- % medidos + erro^2 dos valores estimados
- E1_ARMA = sum((x1ro_s-y1_ARMA2).^2)
- % 3.12) Gráfico do modelo ARMA para 2*tempo
- % estimação de um modelo ARMA
- %tr1_ARMA = polyval(p, tt);
- %tr1_ARMA = 10*sin(tt/1)+20;
- tr1_ARMA = x1;
- st1_ARMA = repmat(y1_ARMA2, 2, 1);
- figure(13)
- plot(t, x1ro, '-+', tt, tr1_ARMA+st1_ARMA, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Modelo ARMA');
- legend('Série temporal', 'Previsão da série para o dobro da duração temporal', 'Location', 'NorthEast');
- % 3.13) Modelo ARIMA da série regularizada
- % definir p, D, q, usando arima e estimate
- p1_ARIMA = 2;
- d1_ARIMA = 1;
- q1_ARIMA = 5;
- Md1 = arima(p1_ARIMA, d1_ARIMA, q1_ARIMA);
- EstMd1 = estimate(Md1, x1ro(1:N), 'Y0', x1ro(1:p1_ARIMA+1));
- y1_ARIMA = simulate(EstMd1, N);
- % 3.14) Teste gráfico usando simulate
- figure(14);
- plot(t, x1ro_s, '-+', t, y1_ARIMA, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Modelo ARIMA');
- legend('Série temporal', 'ARMA calculado com simulate', 'Location', 'NorthEast');
- % 3.15) Validar o modelo comparando a série regularizada com o resultado da
- % combinação aditiva da estimação da componente sazonal + componente da
- % tendencia paramétrica da série. NOTA: mértica = erro^2 dos valores
- % medidos + erro^2 dos valores estimados
- E1_ARIMA = sum((x1ro_s-y1_ARIMA).^2)
- % 3.16) Gráfico do modelo ARIMA para 2*tempo
- %tr1_ARIMA = polyval(p, tt);
- %tr1_ARIMA = 10*sin(tt/1)+20;
- tr1_ARIMA = x1;
- st1_ARIMA = repmat(y1_ARIMA, 2, 1);
- figure(15)
- plot(t, x1ro, '-+', tt, tr1_ARIMA+st1_ARIMA, '-o')
- ylabel('Temperatua(ºC)');
- xlabel('Dia do ano');
- title('Modelo ARIMA');
- legend('Série temporal', 'Previsão da série para o dobro da duração temporal', 'Location', 'NorthEast');
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement