Advertisement
Guest User

Untitled

a guest
Apr 25th, 2018
57
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 13.54 KB | None | 0 0
  1. % 1)
  2. % 1.1) Leitura e representação dos resultados lidos
  3. data = importdata('dataset_ATD_PL7.csv');
  4. leituras = data.data;
  5. datas = data.textdata(2:366);
  6. N = length(leituras);
  7. t = 0:N-1;
  8.  
  9. figure(1);
  10. plot(t, leituras, '-+');
  11. ylabel('Temperatua(ºC)');
  12. xlabel('Dia do ano');
  13. legend('Leituras ao longo do ano','Location', 'NorthEast');
  14.  
  15.  
  16. % 1.2) Verificação da existência de NaN através de extrapolação
  17. x1 = leituras;
  18.  
  19. if any(isnan(x1))
  20. ind = find(isnan(x1));
  21. x1r = x1;
  22. for k=1: length(ind)
  23. tt = t(ind(k)-4 : ind(k)-1);
  24. xx = x1r (ind(k)-4 : ind(k)-1);
  25. x1r(ind(k)) = interp1(tt, xx, t(ind(k)), 'pchip', 'extrap');
  26. end
  27. end
  28.  
  29. figure(2);
  30. plot(t, x1, '-+', t, x1r, '-o');
  31. ylabel('Temperatua(ºC)');
  32. xlabel('Dia do ano');
  33. legend('com NaN', 'sem NaN', 'Location', 'NorthEast');
  34.  
  35. % 1.3) e 1.4) Verificação da existência de outliers através de uma janela deslizante
  36. % valores da média e desvio padrão são calculados a cada novo mês
  37. % substituição dos outliers por valores adequados e sua respetiva
  38. % representação
  39.  
  40. % variáveis auxliares para a janela deslizante
  41. tam = N;
  42. passo=30;
  43. conta=1;
  44. mes=1;
  45.  
  46. % inicializar vetor a 0
  47. s = size(x1r);
  48. x1ro = zeros(s);
  49.  
  50. for i = 1: passo: N
  51. % caso em que o mês tem 30 dias, até ao mês 12
  52. if (mes <= 12)
  53.  
  54. new_x1r = x1r(i:i+(passo-1));
  55. x1ro(i:i+(passo-1)) = new_x1r;
  56.  
  57. % cálculo das médias do mês
  58. mu1 = mean(new_x1r);
  59. sigma1 = std(new_x1r);
  60.  
  61. % verifica 1 a 1 cada dia do mes
  62. for j = 1: 1: 30
  63. if (abs(new_x1r(j) - mu1) > 3*sigma1)
  64.  
  65. disp('outlier')
  66. if x1ro(i+j-1) > mu1
  67. disp(x1ro(i+j-1))
  68. x1ro(i+j-1) = mu1 + 2.5*sigma1;
  69. disp(x1ro(i+j-1))
  70. else
  71. disp(x1ro(i+j-1))
  72. x1ro(i+j-1) = mu1 - 2.5*sigma1;
  73. disp(x1ro(i+j-1))
  74. end
  75. end
  76. end
  77. mes = mes + 1;
  78. % caso do resto dos 5 dias que sobram para o total dos 365 dias
  79. else
  80.  
  81. new_x1r = x1r(i:i+(5-1));
  82. x1ro(i:i+(5-1)) = new_x1r;
  83.  
  84. % cálculo das médias do mês
  85. mu1 = mean(new_x1r);
  86. sigma1 = std(new_x1r);
  87.  
  88. % verifica 1 a 1 cada dia do mes
  89. for j = 1: 1: 5
  90. if (abs(new_x1r(j) - mu1) > 3*sigma1)
  91.  
  92. disp('outlier')
  93. if x1ro(i+j-1) > mu1
  94. disp(x1ro(i+j-1))
  95. x1ro(i+j-1) = mu1 + 2.5*sigma1;
  96. disp(x1ro(i+j-1))
  97. else
  98. disp(x1ro(i+j-1))
  99. x1ro(i+j-1) = mu1 - 2.5*sigma1;
  100. disp(x1ro(i+j-1))
  101. end
  102. end
  103. end
  104. mes = mes + 1;
  105. end
  106. end
  107.  
  108. figure(3);
  109. plot(t, x1r, '-+', t, x1ro, '-o');
  110. ylabel('Temperatua(ºC)');
  111. xlabel('Dia do ano');
  112. legend('com outliers', 'outliers corrigidos', 'Location', 'NorthEast');
  113.  
  114. %% 2
  115. N=length(x1ro);
  116. t=0:N-1;
  117. % 2.1) estimar serie temporal sem a componente de tendencia
  118. % aproximações polinomiais de grau 0 e 1, através da função detrend
  119. % 2.2) obter a componente
  120. x1ro_t0 = detrend(x1ro, 'constant');
  121. x1ro_t1 = detrend(x1ro, 'linear');
  122.  
  123. % 2.3) Representar graficamente a série temporal regularizada, com
  124. % componente e sem componente de tendência
  125. % série temporal com tendência de grau 0
  126. figure(4)
  127. subplot(211)
  128. plot(t, x1ro, '-+', t, x1ro - x1ro_t0, '-*');
  129. ylabel('Temperatua(ºC)');
  130. xlabel('Dia do ano');
  131. title('Série temporal com componente de tendência de grau 0');
  132. legend('Série temporal', 'Componente de tendência', 'Location', 'NorthEast');
  133. subplot(212)
  134. plot(t, x1ro_t0, '-o');
  135. ylabel('Temperatua(ºC)');
  136. xlabel('Dia do ano');
  137. title('Série temporal sem componente de tendência de grau 0');
  138. legend('Sem componente de tendência', 'Location', 'NorthEast');
  139.  
  140. % série temporal com tendência de grau 1
  141. figure(5)
  142. subplot(211)
  143. plot(t, x1ro, '-+', t, x1ro - x1ro_t1, '-*');
  144. ylabel('Temperatua(ºC)');
  145. xlabel('Dia do ano');
  146. title('Série temporal com componente de tendência de grau 1');
  147. legend('Série temporal', 'Componente de tendência', 'Location', 'NorthEast');
  148. subplot(212)
  149. plot(t, x1ro_t1, '-o');
  150. ylabel('Temperatua(ºC)');
  151. xlabel('Dia do ano');
  152. title('Série temporal sem componente de tendência de grau 1');
  153. legend('Sem componente de tendência', 'Location', 'NorthEast');
  154.  
  155. % série com tendência de grau 5, caso polinomial
  156. % 2.4) obter a componente
  157. p = polyfit(t', x1ro, 5);
  158. tr1 = polyval(p, t');
  159.  
  160. % 2.5) Representar graficamente a série temporal com e sem componente de
  161. % tendência
  162. figure(6)
  163. subplot(211)
  164. plot(t, x1ro, '-+', t, tr1, '-*');
  165. ylabel('Temperatua(ºC)');
  166. xlabel('Dia do ano');
  167. title('Série temporal com componente de tendência de grau 5');
  168. legend('Série temporal', 'Com componente de tendência', 'Location', 'NorthEast');
  169. subplot(212)
  170. plot(t, x1ro-tr1, '-o');
  171. ylabel('Temperatua(ºC)');
  172. xlabel('Dia do ano');
  173. title('Série temporal sem componente de tendência de grau 5');
  174. legend('Sem componente de tendência', 'Location', 'NorthEast');
  175.  
  176. %{
  177. % SERIE SINUSOIDAL
  178. t_f = t.*24*3600;
  179. [Cm, tetam] = SerieFourier(t_f, x1ro, 360*24*60*60, 3);
  180. tr1 = 19.8784445884705*sin(2.23762861251355*10^(-17)*t);
  181. mt =[0 1];
  182.  
  183. for k=1: length(mt)
  184. x1 = zeros(size(x));
  185. for m=0:mt(k)
  186. x1 = x1+Cm(m+1)*cos(m*2*pi/T0*t+tetam(m+1));
  187. end
  188. plot(t, x1, '-.b')
  189. end
  190. %}
  191.  
  192. T0 = 365*24*60*60;
  193. m_max = 50;
  194. t_f = t.*24*3600;
  195. [Cm,tetam] = SerieFourier(t_f',x1ro,T0,m_max);
  196. mt =0:50;
  197.  
  198. for k=1:length(mt)
  199. x1 = zeros(size(t));
  200. for m=0:mt(k)
  201. x1 = x1+ Cm(m+1)*cos(m*2*pi/T0*t_f+tetam(m+1));
  202. end
  203. end
  204.  
  205. figure(40)
  206. plot(t_f,x1);
  207. hold on;
  208. tr1 = x1;
  209. %tr1 = 10*sin(t/1)+20;
  210. tr1 = tr1';
  211. figure(20)
  212. subplot(211)
  213. plot(t,tr1,'-',t,x1ro,'-*')
  214. ylabel('Temperatua(ºC)');
  215. xlabel('Dia do ano');
  216. title('Série temporal com componente de tendência sinosuidal');
  217. legend('Série temporal', 'Componente de tendência', 'Location', 'NorthEast');
  218. subplot(212)
  219. plot(t, x1ro-tr1, '-o');
  220. ylabel('Temperatua(ºC)');
  221. xlabel('Dia do ano');
  222. title('Série temporal sem componente de tendência sinosuidal');
  223. legend('Sem componente de tendência', 'Location', 'NorthEast');
  224.  
  225.  
  226. % 2.6) % Estimar a componente da sazonalidade MENSAL (30 dias)
  227. x1ro_t2 = x1ro-tr1;
  228. ho = (1:30)';
  229.  
  230. s = size(ho)*13;
  231. s(2) = 1;
  232. x1ro_s = zeros(s);
  233.  
  234. new_x1ro_t2 = zeros(s);
  235. new_x1ro_t2(1:365) = x1ro_t2(1:365);
  236. new_x1ro_t2(366:end) = x1ro_t2(1:25);
  237. % cálculo da componente da sazonalidade
  238. dv = repmat(ho, 13, 1);
  239. sx = dummyvar(dv);
  240. xtemp = sx\new_x1ro_t2;
  241. x1ro_s = sx * xtemp;
  242.  
  243. x1ro_s = x1ro_s(1:365);
  244.  
  245. % 2.7) Obter série sem a componente da sazonalidade
  246. % 2.8) Representar a série com e sem a componente da sazonalidade
  247. figure(7)
  248. subplot(211)
  249. plot(t, x1ro, '-+', t, x1ro-x1ro_s, '-*', t, tr1, '-o')
  250. ylabel('Temperatua(ºC)');
  251. xlabel('Dia do ano');
  252. title('Série temporal com componente da sazonalidade');
  253. legend('Série temporal', 'Componente da sazonalidade', 'Location', 'NorthEast');
  254. subplot(212)
  255. plot(t, x1ro_s, '-o')
  256. ylabel('Temperatua(ºC)');
  257. xlabel('Dia do ano');
  258. title('Série temporal sem componente da sazonalidade');
  259. legend('Sem da sazonalidade', 'Location', 'NorthEast');
  260.  
  261. % 2.9) obter série temporal com e sem componente irregular
  262. irr_x1ro = x1ro - tr1 - x1ro_s;
  263.  
  264. % 2.10) Representação da serie temporal com e sem componente irregular
  265. figure(8)
  266. subplot(211)
  267. plot(t, x1ro, '-+', t, x1ro-irr_x1ro, '-*')
  268. ylabel('Temperatua(ºC)');
  269. xlabel('Dia do ano');
  270. title('Série temporal com componente irregular');
  271. legend('Com componente irregular', 'Location', 'NorthEast');
  272. subplot(212)
  273. plot(t, irr_x1ro, '-o')
  274. ylabel('Temperatua(ºC)');
  275. xlabel('Dia do ano');
  276. title('Série temporal sem componente irregular');
  277. legend('Sem componente irregular', 'Location', 'NorthEast');
  278.  
  279. %% 3)
  280.  
  281. % 3.1) Obter componentes da série temporal resultantes da decomposição da
  282. % série temporal regularizada
  283. % 3.2) Verificar estacionaridade da série temporal através do adftest
  284. N = length(x1ro);
  285. t = (0:N-1)';
  286. tt = (0:2*N-1)';
  287. adftest(x1ro) % nao é estacionaria = 0
  288. adftest(x1ro_s) % sazonal deu 1- provavelmente é estacionaria
  289. % sazonal, as carateristicas repetem-se igualmente de 24 em 24 horas
  290.  
  291. % componente sazonal da série
  292. y1 = x1ro_s(1:30);
  293.  
  294. % 3.3) Representar graficamente FAC e FACP, através do autocorr e parcorr
  295. figure(9);
  296. subplot(211)
  297. autocorr(y1)
  298. subplot(212)
  299. parcorr(y1)
  300.  
  301. % 3.4) Criar objeto iddata, com Ts=1
  302. id_y1 = iddata(y1, [], 1, 'TimeUnit', 'days');
  303.  
  304. % 3.5) Modelo AR para a componente sazonal, resultado da FACP
  305. % definir na, usando arOptions, ar e polydata
  306. opt1_AR = arOptions('Approach', 'ls');
  307. na1_AR = 10; % repetir o processo ate encontrar numero aconselhavel, mais aceitavel
  308. model1_AR = ar(id_y1, na1_AR, opt1_AR)
  309. pcoef1_AR = polydata(model1_AR);
  310.  
  311. y1_AR = y1(1:na1_AR);
  312.  
  313. for k=na1_AR+1:30 % do ...(na1_AR+1) ao 30
  314. y1_AR(k) = sum(-pcoef1_AR(2:end)'.*flip(y1_AR(k-na1_AR: k-1)));
  315. end
  316.  
  317. % 3.6) Teste gráfico usando parametros do modelo da componente sazonal e forecast
  318. y1_AR2 = repmat(y1_AR, 13, 1);
  319. y1_ARf = forecast(model1_AR, y1(1: na1_AR), 30-na1_AR);
  320. y1_ARf = repmat([y1(1:na1_AR); y1_ARf], 13, 1);
  321.  
  322. % encurtar matrizes para 365 dias
  323. y1_AR2 = y1_AR2(1:365);
  324. y1_ARf = y1_ARf(1:365);
  325.  
  326. figure(10)
  327. plot(t, x1ro_s, '-+', t, y1_AR2, '-o', t, y1_ARf, '-*');
  328. ylabel('Temperatua(ºC)');
  329. xlabel('Dia do ano');
  330. title('Modelo AR');
  331. legend('Série temporal', 'AR calculado com for', 'AR calculado com forecast', 'Location', 'NorthEast');
  332.  
  333. % 3.7) Validar o modelo comparando a série regularizada com o resultado da
  334. % combinação aditiva da estimação da componente sazonal + componente da
  335. % tendencia paramétrica da série. NOTA: mértica = erro^2 dos valores
  336. % medidos + erro^2 dos valores estimados
  337. % calcular erros (métrica)
  338. E1_AR = sum((x1ro_s-y1_AR2).^2);
  339.  
  340. % 3.8) Gráfico do modelo AR para 2*tempo
  341. % aumentar para 2*365 dias
  342. %tr1_AR = polyval(p, tt);
  343. x1 = repmat(x1',2,1);
  344. tr1_AR = x1;
  345. %tr1_AR = 10*sin(tt/1)+20;
  346. st1_AR = repmat(y1_AR2, 2, 1);
  347.  
  348. figure(11)
  349. plot(t, x1ro, '-+', tt, tr1_AR+st1_AR, '-o')
  350. ylabel('Temperatua(ºC)');
  351. xlabel('Dia do ano');
  352. title('Modelo AR');
  353. legend('Série temporal', 'Previsão da série para o dobro da duração temporal', 'Location', 'NorthEast');
  354.  
  355. % 3.9) Modelo ARMA para a componente sazonal, resultado da FAC
  356. % definir na, nc, usando armaxOptions, armax e polydata
  357. opt1_ARMAX = armaxOptions('SearchMethod', 'auto')
  358. ma1_ARMA = 5;
  359. mc1_ARMA = 1;
  360. model1_ARMA = armax(id_y1, [ma1_ARMA mc1_ARMA], opt1_ARMAX);
  361. [pa1_ARMA, pb1_ARMA, pc1_ARMA] = polydata(model1_ARMA);
  362.  
  363. % ruido branco
  364. e = randn(30, 1);
  365. y1_ARMA = y1(1:ma1_ARMA);
  366. for k = ma1_ARMA + 1:30
  367. 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)));
  368. end
  369.  
  370. % 3.10) Teste gráfico usando parametros do modelo da componente sazonal e forecast
  371. % feito o anterior automaticamente com a funcao forcast
  372. y1_ARMA2 = repmat(y1_ARMA, 13, 1);
  373. y1_ARMAf = forecast(model1_ARMA, y1(1:ma1_ARMA), 30-ma1_ARMA);
  374. y1_ARMAf2 = repmat([y1_ARMA(1:ma1_ARMA); y1_ARMAf], 13, 1);
  375.  
  376. % encurtar matrizes para 365 dias
  377. y1_ARMA2 = y1_ARMA2(1:365);
  378. y1_ARMAf2 = y1_ARMAf2(1:365);
  379.  
  380. figure(12);
  381. plot(t, x1ro_s, '-+', t, y1_ARMA2, '-o', t, y1_ARMAf2, '-*');
  382. ylabel('Temperatua(ºC)');
  383. xlabel('Dia do ano');
  384. title('Modelo ARMA');
  385. legend('Série temporal', 'ARMA calculado com for', 'ARMA calculado com forecast', 'Location', 'NorthEast');
  386.  
  387. % 3.11) Validar o modelo comparando a série regularizada com o resultado da
  388. % combinação aditiva da estimação da componente sazonal + componente da
  389. % tendencia paramétrica da série. NOTA: mértica = erro^2 dos valores
  390. % medidos + erro^2 dos valores estimados
  391. E1_ARMA = sum((x1ro_s-y1_ARMA2).^2)
  392.  
  393. % 3.12) Gráfico do modelo ARMA para 2*tempo
  394. % estimação de um modelo ARMA
  395. %tr1_ARMA = polyval(p, tt);
  396. %tr1_ARMA = 10*sin(tt/1)+20;
  397.  
  398. tr1_ARMA = x1;
  399. st1_ARMA = repmat(y1_ARMA2, 2, 1);
  400.  
  401. figure(13)
  402. plot(t, x1ro, '-+', tt, tr1_ARMA+st1_ARMA, '-o')
  403. ylabel('Temperatua(ºC)');
  404. xlabel('Dia do ano');
  405. title('Modelo ARMA');
  406. legend('Série temporal', 'Previsão da série para o dobro da duração temporal', 'Location', 'NorthEast');
  407.  
  408. % 3.13) Modelo ARIMA da série regularizada
  409. % definir p, D, q, usando arima e estimate
  410. p1_ARIMA = 2;
  411. d1_ARIMA = 1;
  412. q1_ARIMA = 5;
  413.  
  414. Md1 = arima(p1_ARIMA, d1_ARIMA, q1_ARIMA);
  415. EstMd1 = estimate(Md1, x1ro(1:N), 'Y0', x1ro(1:p1_ARIMA+1));
  416. y1_ARIMA = simulate(EstMd1, N);
  417.  
  418. % 3.14) Teste gráfico usando simulate
  419. figure(14);
  420. plot(t, x1ro_s, '-+', t, y1_ARIMA, '-o')
  421. ylabel('Temperatua(ºC)');
  422. xlabel('Dia do ano');
  423. title('Modelo ARIMA');
  424. legend('Série temporal', 'ARMA calculado com simulate', 'Location', 'NorthEast');
  425.  
  426. % 3.15) Validar o modelo comparando a série regularizada com o resultado da
  427. % combinação aditiva da estimação da componente sazonal + componente da
  428. % tendencia paramétrica da série. NOTA: mértica = erro^2 dos valores
  429. % medidos + erro^2 dos valores estimados
  430. E1_ARIMA = sum((x1ro_s-y1_ARIMA).^2)
  431.  
  432. % 3.16) Gráfico do modelo ARIMA para 2*tempo
  433. %tr1_ARIMA = polyval(p, tt);
  434. %tr1_ARIMA = 10*sin(tt/1)+20;
  435.  
  436. tr1_ARIMA = x1;
  437. st1_ARIMA = repmat(y1_ARIMA, 2, 1);
  438.  
  439. figure(15)
  440. plot(t, x1ro, '-+', tt, tr1_ARIMA+st1_ARIMA, '-o')
  441. ylabel('Temperatua(ºC)');
  442. xlabel('Dia do ano');
  443. title('Modelo ARIMA');
  444. 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