################################################################################ # Livro "ANÁLISE DE SÉRIES TEMPORAIS" - Morettin e Toloi (2004) ################################################################################ --- & --- # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # INTRODUÇÃO ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # PACOTES require(lmtest) require(tseries) require(forecast) ################################################################################ # Explorando a Série ICV ################################################################################ dados5 <- read.table('D-BANESPA.txt');dados5 ICV <- read.table('D-BANESPA.txt', head=F) dim(icv) icv <- ICV[,1] summary(icv) ts.plot(icv) icv<-ICV[,1] plot(icv) plot(icv, type="l") plot(icv, type="l", lwd=2) plot(icv, type="l", lwd=1.5, col="blue") plot(icv, type="l", lwd=2, col="blue", xlab="Anos") plot(icv, type="l", lwd=2, col="blue", xlab="Anos", ylab="ICV", ylim = c(10,800)) ################################################################################ # Testing Normality ################################################################################ dev.off() x11() hist(icv, prob=T, col="red") lines(density(icv), lwd=2) ################################################################################ # Fitting a regression model ################################################################################ t<-seq(1970,1979,length=length(icv)) t t2<-t^2 reg<-lm(icv~t+t2) reg names(reg) reg$fit anova(reg) summary(reg) plot(icv, type="l") lines(reg$fit,col=2,lwd=2) resicv<-resid(reg) plot(resicv) ################################################################################ # Checking the variability ################################################################################ icv1 <- (icv[1:8]) icv1 m1 <- mean(icv1) d1 <- var(icv1)^.5 w1 <- (max(icv1)-min(icv1)) icv2 <- (icv[9:16]) icv2 m2 <- mean(icv2) d2 <- var(icv2)^.5 w2 <- (max(icv2)-min(icv2)) icv3 <- (icv[17:24]) icv3 m3 <- mean(icv3) d3 <- var(icv3)^.5 w3 <- (max(icv3)-min(icv3)) icv4 <- (icv[25:32]) icv4 m4 <- mean(icv4) d4 <- var(icv4)^.5 w4 <- (max(icv4)-min(icv4)) icv5 <- (icv[33:40]) icv5 m5 <- mean(icv5) d5 <- var(icv5)^.5 w5 <- (max(icv5)-min(icv5)) icv6 <- (icv[41:48]) icv6 m6 <- mean(icv6) d6 <- var(icv6)^.5 w6 <- (max(icv6)-min(icv6)) icv7 <- (icv[49:56]) icv7 m7 <- mean(icv7) d7 <- var(icv7)^.5 w7 <- (max(icv7)-min(icv7)) icv8 <- (icv[57:64]) icv8 m8 <- mean(icv8) d8 <- var(icv8)^.5 w8 <- (max(icv8)-min(icv8)) icv9 <- (icv[65:72]) icv9 m9 <- mean(icv9) d9 <- var(icv9)^.5 w9 <- (max(icv9)-min(icv9)) icv91 <- (icv[73:80]) icv91 m91 <- mean(icv91) d91 <- var(icv91)^.5 w91 <- (max(icv91)-min(icv91)) icv92 <- (icv[81:88]) icv92 m92 <- mean(icv92) d92 <- var(icv92)^.5 w92 <- (max(icv92)-min(icv92)) icv93 <- (icv[89:96]) icv93 m93 <- mean(icv93) d93 <- var(icv93)^.5 w93 <- (max(icv93)-min(icv93)) icv94 <- (icv[97:104]) icv94 m94 <- mean(icv94) d94 <- var(icv94)^.5 w94 <- (max(icv94)-min(icv94)) icv95 <- (icv[105:114]) icv95 m95 <- mean(icv95) d95 <- var(icv95)^.5 w95 <- (max(icv95)-min(icv95)) ZbarDP <- cbind(rbind(m1, m2, m3, m4, m5, m6, m7, m8, m9, m91, m92, m93, m94, m95),rbind(d1, d2, d3, d4, d5, d6, d7, d8, d9, d91, d92, d93, d94, d95)) ZbarDP ZbarW <- cbind(rbind(m1, m2, m3, m4, m5, m6, m7, m8, m9, m91, m92, m93, m94, m95),rbind(w1, w2, w3, w4, w5, w6, w7, w8, w9, w91, w92, w93, w94, w95)) ZbarW par(mfrow = c(1,2)) plot(ZbarDP) plot(ZbarW) ################################################################################ # Checking the variability - moving average ################################################################################ dev.off() x11() j <- floor(114/12);j n <- 114 abscissa <- NULL ordenada <- NULL for (i in 1: (n-j)){ abscissa[i] <- mean(icv[i:(i+j)]) ordenada[i] <- sd(icv[i:(i+j)]) } plot(x=abscissa,y=ordenada, xlab='Média', ylab= 'Desvio-padrão', main='Verificação da Variabilidade') ################################################################################ # Checking the variability - grouping (for large samples) ################################################################################ n <- 114 # total number of observations g <- 12 # number of observations in each group g1 <- g-1 j <- floor(n/g);j # number of groups SEQ <- seq (1, n, by=g) abscissa1 <- NULL ordenada1 <- NULL for (i in SEQ){ abscissa1[i] <- mean(icv[i:(i+g1)]) ordenada1[i] <- sd(icv[i:(i+g1)]) } abscissa1 ordenada1 abscissa11 <- abscissa1[1:97] ordenada11 <- ordenada1[1:97] plot(x=abscissa11,y=ordenada11, xlab='Média', ylab= 'Desvio-padrão', main='Verificação da Variabilidade') # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # MODELOS ARIMA (p,d,q) ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ------------------------------------------------------------------------------ # Modelos Auto-regressivos - AR(p) # ------------------------------------------------------------------------------ # AR(1) com phi = 0.6 ar.sim <- arima.sim(n=10000, list(ar=c(0.6))) ar.sim plot(ar.sim, main="processo AR(1) simulado com phi=0.6", type="l") par(mfrow=c(1,2)) acf.ar <- acf(ar.sim) pacf(ar.sim) acf.ar$acf[2] # AR(1) com phi = -0.6 simulacao <- arima.sim(n=10000, list(ar=c(-0.6))) acf(simulacao) pacf(simulacao) acf.ar$acf[2] # AR(1) com phi = 0.1 simulacao <- arima.sim(n=10000, list(ar=c(0.1))) acf(simulacao) pacf(simulacao) acf.ar$acf[2] # AR(1) com phi = 0.01 simulacao <- arima.sim(n=10000, list(ar=c(0.01))) ar.sim acf(simulacao) pacf(simulacao) acf.ar$acf[2] # AR(1) com phi = 0.9 simulacao <- arima.sim(n=10000, list(ar=c(0.9))) acf(simulacao) pacf(simulacao) acf.ar$acf[2] # AR(1) com phi = 0.99 simulacao <- arima.sim(n=10000, list(ar=c(0.99))) acf(simulacao) pacf(simulacao) acf.ar$acf[2] # AR(2) com phi_1 = 0.8 e phi_2 = -0.5 ar.sim <- arima.sim(n=10000, list(ar=c(0.3,-0.5))) acf(simulacao) pacf(simulacao) # AR(2) com phi_1 = -0.8 e phi_2 = 0.5 simulacao <- arima.sim(n=10000, list(ar=c(-0.3,0.5))) acf(simulacao) pacf(simulacao) # AR(2) com phi_1 = -0.8 e phi_2 = 0.5 # verificar condicoes de estacionariedade ar.sim <- arima.sim(n=10000, list(ar=c(-0.8,0.5))) acf(simulacao) pacf(simulacao) dev.off() x11() # ------------------------------------------------------------------------------ # Modelos de Médias Móveis - MA(q) # ------------------------------------------------------------------------------ ?arima.sim # MA(1) com theta = 0.8 # Obs: o R considera o modelo com -theta ma.sim <- arima.sim(n=10000, list(ma=c(-0.8))) ma.sim plot(ma.sim, main="processo MA(1) simulado com theta=-0.8", type="l") par(mfrow=c(1,2)) acf.ma <- acf(ma.sim) pacf(ma.sim) acf.ma$acf[2] # lembrando que ro_j = -theta/(1+theta^2) ro_j <- -0.8/(1+0.8^2) ro_j # MA(1) com theta = 0.1 ma.sim <- arima.sim(n=10000, list(ma=c(-0.1))) ma.sim acf.ma <- acf(ma.sim) pacf(ma.sim) # MA(1) com theta = 0.9 ma.sim <- arima.sim(n=10000, list(ma=c(-0.9))) ma.sim acf.ma <- acf(ma.sim) pacf(ma.sim) # MA(1) com theta = 0.99 ma.sim <- arima.sim(n=10000, list(ma=c(-0.99))) ma.sim acf.ma <- acf(ma.sim) pacf(ma.sim) # MA(2) com theta_1 = 0.5 e theta_2 = -0.3 ma.sim <- arima.sim(n=10000, list(ma=c(-0.5,0.3))) ma.sim acf.ma <- acf(ma.sim) pacf(ma.sim) # MA(2) com theta_1 = -0.5 e theta_2 = 0.3 ma.sim <- arima.sim(n=10000, list(ma=c(0.5,-0.3))) ma.sim acf.ma <- acf(ma.sim) pacf(ma.sim) # MA(2) com theta_1 = -0.5 e theta_2 = 0.8 # verificar condicoes de invertibilidade # Cuidado!! O R não acusa processos MA(q) NÃO invertíveis. ma.sim <- arima.sim(n=10000, list(ma=c(0.5,-0.8))) ma.sim acf.ma <- acf(ma.sim) pacf(ma.sim) # ------------------------------------------------------------------------------ # Modelos Auto-regressivos e de Médias Móveis - ARMA(p,q) # ------------------------------------------------------------------------------ # ARMA(1,1) com phi = 0.8 e theta = 0.3 arma.sim <- arima.sim(n=10000, list(ar=c(0.8), ma=c(-0.3))) arma.sim ts.plot(arma.sim) acf.arma <- acf(arma.sim) pacf(arma.sim) acf.arma # ARMA(2,1) com phi_1 = 0.7, phi_2 = 0.2 e theta = 0.3 arma.sim <- arima.sim(n=10000, list(ar=c(0.7,0.2), ma=c(0.3))) arma.sim ts.plot(arma.sim) acf.arma <- acf(arma.sim) pacf(arma.sim) # ------------------------------------------------------------------------------ # Modelos Auto-regressivos Integrados e de Médias Móveis - ARIMA(p,d,q) # ------------------------------------------------------------------------------ # ARIMA(1,1,0) com phi = 0.8 ts.sim <- arima.sim(list(order = c(1,1,0), ar = 0.4), n = 1000) ts.plot(ts.sim) acf.arima <- acf(ts.sim) pacf(ts.sim) acf(diff(ts.sim)) pacf(diff(ts.sim)) # ARIMA(0,1,1) com theta = 0.3 ts.sim <- arima.sim(list(order = c(0,1,1), ma = -0.3), n = 1000) ts.plot(ts.sim) acf.arima <- acf(ts.sim) pacf(diff(ts.sim)) acf(diff(ts.sim)) pacf(diff(ts.sim)) # ARIMA(1,1,1) com phi = 0.8 e theta = 0.3 ts.sim <- arima.sim(list(order = c(1,1,1), ar = 0.8, ma = -0.3), n = 1000) ts.plot(ts.sim) acf.arima <- acf(ts.sim) pacf(diff(ts.sim)) acf(diff(ts.sim)) pacf(diff(ts.sim)) # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # Identificação de Modelos ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # Identificando série simuladas # Obs: O R calcula o IC baseado em uma série não correlacionada, por isso deve # ser interpretada com cuidado!! # Sugiro usar a opção ci.type = "ma" para um processo em que as autocorrelações # são nulas para v > q (MA(q)). ar.sim <- arima.sim(n=10, list(ar=c(0.3))) par(mfrow=c(1,2)) acf.ar <- acf(ar.sim) acf.ar <- acf(ar.sim, ci.type="ma") pacf(ar.sim) plot(acf.ar,ci.type="ma",ci=0.95) plot(acf.ar,ci.type="ma",ci=0.90) plot(acf.ar,ci.type="ma",ci=0.99) par(mfrow=c(1,2)) simuladas <- arima.sim(n=10000, list(ar=c(0.3))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(-0.8))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(-0.2, 0.5))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(0.2, -0.2, 0.5))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(0.2, -0.2, -0.5, 0.1))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(0.5, 0.5))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(0.9))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(0.999))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ma=c(-0.2))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ma=c(-0.2, 0,3))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") simuladas <- arima.sim(n=10000, list(ar=c(0.7,0.1), ma=c(-0.2, 0.5))) acf(simuladas, ci.type="ma", xlab = "defasagens", ylab = "FAC") pacf(simuladas, xlab = "defasagens", ylab = "FAC parcial") ################################################################################ # Exemplo 6.3 - série ICV ################################################################################ icv plot(icv, main="ICV de jan-70 a jun-79", lwd = "2", col = "51", xlab = "anos", ylim = c(10,85)) summary(icv) length(icv) # ------------------------------------------------------------------------------ # Verificação da Variabilidade para n pequeno(médias móveis) # ------------------------------------------------------------------------------ n <- 114 j <- floor(n/13);j abscissa <- NULL ordenada <- NULL for (i in 1: (n-j)){ abscissa[i] <- mean(icv[i:(i+j)]) ordenada[i] <- sd(icv[i:(i+j)]) } plot(abscissa,ordenada, xlab='Média', ylab= 'Desvio-padrão', main='Verificação da Variabilidade') # ------------------------------------------------------------------------------ # Verificação da Variabilidade para n grande (intervalo de faixa) # ------------------------------------------------------------------------------ n <- 114 # número de observações g <- 8 # número de observações em cada grupo n1 <- n-g g1 <- g-1 j <- floor(n/g) # número de grupos SEQ <- seq (1, n1, by=g) abscissa1 <- NULL ordenada1 <- NULL for (i in SEQ){ abscissa1[i] <- mean(icv[i:(i+g1)]) ordenada1[i] <- sd(icv[i:(i+g1)]) } plot(x=abscissa1,y=ordenada1, xlab='Média', ylab= 'Desvio-padrão', main='Verificação da Variabilidade') # ------------------------------------------------------------------------------ # Transformando os dados # ------------------------------------------------------------------------------ logicv <- log(icv) par(mfrow=c(2,1)) ts.plot(icv, type = "l", main="série sem transformação") ts.plot(logicv, type = "l", main="série com transformação") # ------------------------------------------------------------------------------ # Identificando o modelo # ------------------------------------------------------------------------------ par(mfrow=c(3,2)) acf(logicv, main = "log icv") pacf(logicv, main = "log icv") d1logicv <- diff(logicv) acf(d1logicv, main = "1a diferença do log icv") pacf(d1logicv, main = "1a diferença do log icv") d2logicv <- diff(d1logicv) acf(d2logicv, main = "2a diferença do log icv") pacf(d2logicv, main = "2a diferença do log icv") # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(1,1,0) com uma constante # (1 - B)(1 - phi_1*B)Yt = theta_0 + at, com Yt = log(icv) # ------------------------------------------------------------------------------ ################################################################################ # Exemplo 6.4 - série Umidade ################################################################################ UMID <- read.csv("C:\\R\\atmosfera.csv", head=T) UMID umid <- UMID [1:358,3] temp <- UMID [,2] summary(umid) ts.plot(umid, lwd = "2", col = "1", xlab = "dias") length(umid) # ------------------------------------------------------------------------------ # Verificação da Variabilidade para n grande (intervalo de faixa) # ------------------------------------------------------------------------------ n <- 358 # total number of observations g <- 7 # number of observations in each group g1 <- g-1 j <- floor(n/g) # number of groups SEQ <- seq (1, n, by=g) abscissa1 <- NULL ordenada1 <- NULL for (i in SEQ){ abscissa1[i] <- mean(umid[i:(i+g1)]) ordenada1[i] <- sd(umid[i:(i+g1)]) } abscissa1 ordenada1 abscissa11 <- abscissa1[1:358] ordenada11 <- ordenada1[1:358] plot(x=abscissa11,y=ordenada11, xlab='Média', ylab= 'Desvio-padrão', main='Verificação da Variabilidade') # ------------------------------------------------------------------------------ # Aparentemente sem necessidade de transformar os dados # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ # Identificando o modelo # ------------------------------------------------------------------------------ par(mfrow=c(3,2)) acf(umid) pacf(umid) d1umid <- diff(umid) acf(d1umid, main = "1a diferença") pacf(d1umid, main = "1a diferença") d2umid <- diff(d1umid) acf(d2umid, main = "2a diferença") pacf(d2umid, main = "2a diferença") # ------------------------------------------------------------------------------ # Quantas diferenças? # ------------------------------------------------------------------------------ # veja a FAC - 1o lag: ~ -0.5 # veja a FACP - mudança de sinais dos valores da FACP # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(1,0,0) com uma constante # (1 - phi_1*B)Zt = theta_0 + at # ------------------------------------------------------------------------------ ################################################################################ # Exemplo 6.5 - série CO ################################################################################ CO <- read.csv("C:\\R\\poluicao.csv", head=T) CO names(CO) co <- CO [1:358,3];co summary(co) ts.plot(co, lwd = "2", col = "1", xlab = "dias") length(co) # ------------------------------------------------------------------------------ # Verificação da Variabilidade para n grande (intervalo de faixa) # ------------------------------------------------------------------------------ n <- 358 # total number of observations g <- 7 # number of observations in each group g1 <- g-1 j <- floor(n/g) # number of groups SEQ <- seq (1, n, by=g) abscissa1 <- NULL ordenada1 <- NULL for (i in SEQ){ abscissa1[i] <- mean(co[i:(i+g1)]) ordenada1[i] <- sd(co[i:(i+g1)]) } abscissa1 ordenada1 abscissa11 <- abscissa1[1:358] ordenada11 <- ordenada1[1:358] plot(x=abscissa11,y=ordenada11, xlab='Média', ylab= 'Desvio-padrão', main='Verificação da Variabilidade') # ------------------------------------------------------------------------------ # Transformação dos dados # ------------------------------------------------------------------------------ logco <- log(co) ts.plot(logco) # ------------------------------------------------------------------------------ # Identificando o modelo # ------------------------------------------------------------------------------ par(mfrow=c(3,2)) acf(logco) pacf(logco) d1logco <- diff(logco) acf(d1logco, main = "1a diferença") pacf(d1logco, main = "1a diferença") d2logco <- diff(d1logco) acf(d2logco, main = "2a diferença") pacf(d2logco, main = "2a diferença") # ------------------------------------------------------------------------------ # Quantas diferenças? # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(3,0,0) com uma constante, com Yt = log(icv) # ------------------------------------------------------------------------------ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # ESTIMAÇÃO DE MODELOS ARIMA(p,d,q) ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ################################################################################ # Exemplo 7.3 Estimativa dos parâmetros das séries simuladas ################################################################################ # ------------------------------------------------------------------------------ # Processo AR(1) # ------------------------------------------------------------------------------ ar.sim <- arima.sim(n=10000, list(ar=c(0.8))) ar_fit <- arima(ar.sim, order=c(1,0,0), include.mean=F);ar_fit #significância do parâmetro - altamente significativo 0.7936/0.0061 # ------------------------------------------------------------------------------ # Processo MA(1) # ------------------------------------------------------------------------------ ma.sim <- arima.sim(n=10000, list(ma=c(-0.8))) ma_fit <- arima(ma.sim, order=c(0,0,1), include.mean=F);ma_fit #significância do parâmetro - altamente significativo -0.7919/0.0059 coeftest(ma_fit) # ------------------------------------------------------------------------------ # Processo ARMA(1,1) # ------------------------------------------------------------------------------ arma.sim <- arima.sim(n=10000, list(ar=c(0.8), ma=c(-0.3))) arma_fit <- arima(arma.sim, order=c(1,0,1), include.mean=F);arma_fit #significância dos parâmetros - altamente significativos 0.8043/0.009 -0.2954/0.0146 ################################################################################ # Exemplo 7.4 Retornando ao exemplo 6.3 - série ICV ################################################################################ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(1,1,0) com uma constante # (1 - B)(1 - phi_1*B)Yt = theta_0 + at, com Yt = log(icv) # ------------------------------------------------------------------------------ # arima_icv <- arima(logicv, order=c(1,1,0)); arima_icv # A função "arima" não permite incluir uma constante quando d=1! Optar pela # opção "xreg"!! theta0 <- seq(1,length(logicv),1) arima_icv <- arima(logicv,order=c(1,1,0), xreg = theta0);arima_icv # O que o R chama de "intercept" nos resultados é a média do processo e não a # constante (intercepto). coeftest(arima_icv) arima_icv1 <- arima(diff(logicv), order=c(1,0,0), include.mean=T); arima_icv1 # Podemos optar por x = diff(logicv), mas a previsão resultará em valores # diferençados!! coeftest(arima_icv1) # ------------------------------------------------------------------------------ # Cálculo da constante # ------------------------------------------------------------------------------ theta0 <- 0.0014*(1-0.1358); theta0 # ------------------------------------------------------------------------------ # Cálculo do erro padrão # ------------------------------------------------------------------------------ ep <- (arima_icv1$var.coef)^0.5;ep # ------------------------------------------------------------------------------ # Significância dos parâmetros # ------------------------------------------------------------------------------ # phi1 arima_icv$coef[1] ep[1,1] arima_icv$coef[1]/ep[1,1] # intercept arima_icv$coef[2] ep[2,2] arima_icv$coef[2]/ep[2,2] # ------------------------------------------------------------------------------ # Outra forma de estimar modelos ARIMA(p,d,q) # ------------------------------------------------------------------------------ arma_icv <- arma(d1logicv,order=c(1,0), include.intercept=T); arma_icv # nessa função o resultado "intercept" é de fato a constante! coeftest(arma_icv) ################################################################################ # Exemplo 7.5 Retornando ao exemplo 6.4 - série UMIDADE ################################################################################ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(1,0,0) com uma constante # (1 - phi_1*B)Zt = theta_0 + at # ------------------------------------------------------------------------------ arima_umid <- arima(umid, order=c(1,0,0), include.mean=T); arima_umid coeftest(arima_umid) # ------------------------------------------------------------------------------ # Cálculo da constante # ------------------------------------------------------------------------------ theta0 <- arima_umid$coef[2]*(1-arima_umid$coef[1]); theta0 # ------------------------------------------------------------------------------ # Cálculo do erro padrão # ------------------------------------------------------------------------------ ep <- (arima_umid$var.coef)^0.5;ep # ------------------------------------------------------------------------------ # Significância dos parâmetros # ------------------------------------------------------------------------------ # phi1 arima_umid$coef[1]/ep[1,1] # intercept arima_umid$coef[2]/ep[2,2] # ------------------------------------------------------------------------------ # Outra forma de estimar modelos ARIMA(p,d,q) # ------------------------------------------------------------------------------ arma_umid <- arma(umid,order=c(1,0), include.intercept=T); arma_umid ################################################################################ # Exemplo 7.6 Retornando ao exemplo 6.5 - série CO ################################################################################ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(3,0,0) com uma constante # ------------------------------------------------------------------------------ arima_co <- arima(logco, order=c(3,0,0), include.mean=T); arima_co coeftest(arima_co) # ------------------------------------------------------------------------------ # Cálculo da constante # ------------------------------------------------------------------------------ theta0 <- arima_co$coef[4]*(1-arima_co$coef[1]-arima_co$coef[2]-arima_co$coef[3]) theta0 # ------------------------------------------------------------------------------ # Cálculo do erro padrão # ------------------------------------------------------------------------------ ep <- (arima_co$var.coef)^0.5;ep # ------------------------------------------------------------------------------ # Significância dos parâmetros # ------------------------------------------------------------------------------ # phi1 arima_co$coef[1]/ep[1,1] # phi2 arima_co$coef[2]/ep[2,2] # phi3 arima_co$coef[3]/ep[3,3] # intercept arima_co$coef[4]/ep[4,4] # ------------------------------------------------------------------------------ # Outra forma de estimar modelos ARIMA(p,d,q) # ------------------------------------------------------------------------------ arma_co <- arma(logco,order=c(3,0), include.intercept=T); arma_co # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # DIAGNÓSTICO DE MODELOS ARIMA ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ################################################################################ # Exemplo 8.2 Retornando ao exemplo 6.3 - série ICV ################################################################################ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(1,1,0) # (1 - B)(1 - phi_1*B)Yt = at, com Yt = log(icv) # ------------------------------------------------------------------------------ arima_icv <- arima(logicv,order=c(1,1,0));arima_icv #arima_icv1 <- arima(diff(logicv), order=c(1,0,0), include.mean=T); arima_icv1 #arma_icv <- arma(d1logicv,order=c(1,0), include.intercept=T); arma_icv # ------------------------------------------------------------------------------ # Análise dos resíduos: modelo "arima_icv" # ------------------------------------------------------------------------------ tsdiag(arima_icv) par(mfrow=c(2,1)) res_arima_icv <- arima_icv$residuals hist(res_arima_icv,freq=F) qqnorm(res_arima_icv) qqline(res_arima_icv) # ------------------------------------------------------------------------------ # H0: os dados não são normais # ------------------------------------------------------------------------------ shapiro.test(res_arima_icv) # ------------------------------------------------------------------------------ # Análise do acf e pacf # ------------------------------------------------------------------------------ acf(arima_icv$resid, ci.type="ma") pacf(arima_icv$resid) # Verificamos que houve um bom ajuste aos dados # ------------------------------------------------------------------------------ # Teste Box-Pierce # H0 null hypothesis of independence in a given time series # ------------------------------------------------------------------------------ Box.test(arima_icv$residuals,lag=10, type="Box-Pierce") # ------------------------------------------------------------------------------ # Teste Ljung-Box # ------------------------------------------------------------------------------ tsdiag(arima_icv) for (i in 1:10){ valorP = Box.test(arima_icv$residuals,i,type='Ljung-Box')$p.value print(valorP) } # ------------------------------------------------------------------------------ # Cálculo da variância residual # ------------------------------------------------------------------------------ sigma2_hat_res_icv <- var(arima_icv$residuals) sigma2_hat_res_icv ################################################################################ # Exemplo 8.3 Retornando ao exemplo 6.4 - série UMIDADE ################################################################################ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(1,0,0) com uma constante # (1 - phi_1*B)Zt = theta_0 + at # ------------------------------------------------------------------------------ arima_umid <- arima(umid, order=c(1,0,0), include.mean=T); arima_umid #arma_umid <- arma(umid,order=c(1,0), include.intercept=T); arma_umid # ------------------------------------------------------------------------------ # Análise dos resíduos: modelo "arima_umid" # ------------------------------------------------------------------------------ tsdiag(arima_umid) par(mfrow=c(2,1)) res_arima_umid <- arima_umid$residuals hist(res_arima_umid,freq=F) qqnorm(res_arima_umid) qqline(res_arima_umid) # ------------------------------------------------------------------------------ # H0: os dados não são normais # ------------------------------------------------------------------------------ shapiro.test(res_arima_umid) # ------------------------------------------------------------------------------ # Análise do acf e pacf # ------------------------------------------------------------------------------ acf(arima_umid$resid, ci.type="ma") # verificamos que r2 é significante pacf(arima_umid$resid) # verificamos que phi22 é significante # ccf(umid,arima_umid$resid) # Verificamos que o modelo não é adequado! # ------------------------------------------------------------------------------ # Vamos modificar o modelo incluindo um parâmetro AR de 2º ordem. # ------------------------------------------------------------------------------ arima_umid2 <- arima(umid, order=c(2,0,0), include.mean=T); arima_umid2 #arma_umid2 <- arma(umid,order=c(2,0), include.intercept=T); arma_umid2 # ------------------------------------------------------------------------------ # Plotando o acf e pacf residual do novo modelo ajustado # ------------------------------------------------------------------------------ acf(arima_umid2$residuals,na.action = na.pass, ci.type="ma") # verificamos que r15 é significante. pacf(arima_umid2$residuals,na.action = na.pass) # verificamos que phi15,15 é significante. Vamos introduzir um parâmetro phi15 # ------------------------------------------------------------------------------ # Novo Modelo # ------------------------------------------------------------------------------ arma_umid3 <- arma(umid,lag=list(ar=c(1,2,15)),include.intercept=T); arma_umid3 # plotando o acf e pacf residual do novo modelo ajustado par(mfrow=c(1,2)) acf(arma_umid3$residuals,na.action = na.pass, ci.type="ma") pacf(arma_umid3$residuals,na.action = na.pass) # verificamos que o modelo incompleto é adequado! Verificamos que as fac e facp # dos resíduos do modelo não indicam nenhuma quebra de comportamento de RB nos # resíduos. # Outra forma de estimar modelos incompletos - usando a opção "fixed"!! arima_umid <- arima(umid, order=c(15,0,0),fixed=c(NA,NA,0,0,0,0,0,0,0,0,0,0,0,0, NA,NA),include.mean=T);arima_umid # Obs: fixed = (phi, PHI, theta, THETA, mean) - coef. simples em minúsculo e # coef. sazonais em maiúsculo acf(arima_umid$residuals,na.action = na.pass, ci.type="ma") pacf(arima_umid$residuals,na.action = na.pass) # ------------------------------------------------------------------------------ # Teste Box-Pierce # H0 null hypothesis of independence in a given time series # ------------------------------------------------------------------------------ Box.test(arma_umid3$residuals,lag=10, type="Box-Pierce") # ------------------------------------------------------------------------------ # Teste Ljung-Box # ------------------------------------------------------------------------------ tsdiag(arima_umid) for (i in 1:10){ valorP = Box.test(arma_umid3$residuals,i,type='Ljung-Box')$p.value print(valorP) } # ------------------------------------------------------------------------------ # Cálculo da variância residual # ------------------------------------------------------------------------------ sigma2_hat_res_umid <- var(arma_umid3$residuals[16:358]) sigma2_hat_res_umid ################################################################################ # Exemplo 8.4 Retornando ao exemplo 6.5 - série CO ################################################################################ # ------------------------------------------------------------------------------ # UM MODELO PRELIMINAR É ARIMA(3,0,0) com uma constante, com Yt = log(icv) # ------------------------------------------------------------------------------ arma_co <- arma(logco,order=c(3,0), include.intercept=T) summary(arma_co) arima_co <- arima(logco, order=c(3,0,0), include.mean=T);arima_co # ------------------------------------------------------------------------------ # Teste Box-Pierce # H0 null hypothesis of independence in a given time series # ------------------------------------------------------------------------------ Box.test(arima_co$residuals,lag=24, type="Box-Pierce") # ------------------------------------------------------------------------------ # Teste Ljung-Box # ------------------------------------------------------------------------------ for (i in 1:50){ valorP <- Box.test(arima_co$residuals,i,type='Ljung-Box')$p.value print(valorP) } # ------------------------------------------------------------------------------ # Análise dos resíduos: modelo "arima_logco" # ------------------------------------------------------------------------------ tsdiag(arima_co) par(mfrow=c(2,1)) res_arima_co <- arima_co$residuals hist(res_arima_co,freq=F) qqnorm(res_arima_co) qqline(res_arima_co) # ------------------------------------------------------------------------------ # H0: os dados não são normais # ------------------------------------------------------------------------------ shapiro.test(res_arima_co) # ------------------------------------------------------------------------------ # Análise do acf e pacf # ------------------------------------------------------------------------------ acf(arima_co$resid, ci.type="ma") # verificamos que r14,r16, r21,r25 são significantes pacf(arima_co$resid) # verificamos que r14,r16, r21,r25 são significantes # ccf(umid,arima_umid$resid) # Verificamos que o modelo não é adequado! # ------------------------------------------------------------------------------ # Vamos modificar o modelo incluindo parâmetros MA de ordens 14,16,21,25. # ------------------------------------------------------------------------------ # não se esqueça da função "arima" com a opção "fixed" para estimar modelos # incompletos!!!!!!!!!!!!! arma_co2 <- arma(logco,lag=list(ar=c(1,2,3), ma=c(14,16,21,25)), include.intercept=T) summary(arma_co2) # plotando o acf e pacf residual do novo modelo ajustado acf(arma_co2$residuals,na.action = na.pass, ci.type="ma") pacf(arma_co2$residuals,na.action = na.pass) # verificamos que o modelo incompleto é adequado! Verificamos que as fac e facp # dos resíduos do modelo não indicam nenhuma quebra de comportamento de RB nos # resíduos. # ------------------------------------------------------------------------------ # Podemos também optar por parâmetros AR de ordens 14,16,21,25. # ------------------------------------------------------------------------------ arma_co3 <- arma(logco,lag=list(ar=c(1,2,3,14,16,21,25)),include.intercept=T) summary(arma_co3) # plotando o acf e pacf residual do novo modelo ajustado acf(arma_co3$residuals,na.action = na.pass, ci.type="ma") pacf(arma_co3$residuals,na.action = na.pass) # ------------------------------------------------------------------------------ # Observe o critério AIC e note que o modelo com parâmetros AR é preferível ao # modelo com parâmetros MA. # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ # Teste Box-Pierce # H0 null hypothesis of independence in a given time series # ------------------------------------------------------------------------------ Box.test(arma_co3$residuals,lag=12, type="Box-Pierce") Box.test(arma_co3$residuals,lag=24, type="Box-Pierce") # ------------------------------------------------------------------------------ # Teste Ljung-Box # ------------------------------------------------------------------------------ for (i in 1:36){ valorP = Box.test(arma_co3$residuals,i,type='Ljung-Box')$p.value print(valorP) } # ------------------------------------------------------------------------------ # Cálculo da variância residual # ------------------------------------------------------------------------------ sigma2_hat_res_co <- var(arma_co3$residuals[26:358]) sigma2_hat_res_co # ------------------------------------------------------------------------------ # Um pouco mais sobre diagnóstico...MODELOS BASEADOS EM CRITERIOS DE ESCOLHA # DE MODELOS # ------------------------------------------------------------------------------ # ------------------------------------------------------------------------------ # Série ICV # ------------------------------------------------------------------------------ # Modelos ajustado arima_icv <- arima(logicv,order=c(1,1,0), xreg = theta0);arima_icv # Série ICV (ARIMA (1,2,1), com cte) arima_icv11 <- arima(logicv,order=c(1,2,1), xreg = theta0);arima_icv11 # Série ICV (ARIMA (1,0,2), com cte) arima_icv12 <- arima(logicv,order=c(1,0,2), xreg = theta0);arima_icv12 # Série ICV (ARIMA (2,1,1), com cte) arima_icv13 <- arima(logicv,order=c(2,1,1), xreg = theta0);arima_icv13 # Série ICV (ARIMA (0,2,1), sem cte) arima_icv14 <- arima(logicv,order=c(0,2,1));arima_icv14 # Série ICV (ARIMA (1,1,1), sem cte) arima_icv15 <- arima(logicv,order=c(1,1,1));arima_icv15 AICs <- cbind(arima_icv$aic, arima_icv11$aic,arima_icv12$aic,arima_icv13$aic, arima_icv14$aic,arima_icv15$aic);AICs acf(arima_icv11$residuals) pacf(arima_icv11$residuals) # este modelo apresenta bom ajuste e minimiza o critério AIC!! auto <- auto.arima(logicv) auto A função "auto.arima" recomenda o uso deste modelo!! # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 9 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # PREVISÃO COM MODELOS ARIMA ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # O R possui duas funções: "forecast" e "predict" ################################################################################ # Exemplo 9.11 Retornando ao exemplo 6.3 - série ICV ################################################################################ arima_icv <- arima(logicv, order=c(1,1,0)); arima_icv require(forecast) prev_icv <- forecast(arima_icv, h = 12); prev_icv plot(prev_icv) names(prev_icv) prev_icv$mean exp(prev_icv$mean) prev_icv_auto <- forecast(auto, h = 12) prev_icv_auto plot(prev_icv_auto) ################################################################################ # Exemplo 9.11 Retornando ao exemplo 6.3 - série UMID ################################################################################ #require(farma) #arma_umid1 <- armaFit(umid, fixed = c(0.1,0.4), data = umid,method = c("mle"), #include.mean=TRUE) #arma_umid1 #prev_umid1 <- predict(arma_umid1, n.ahead = 12); prev_umid1 #plot(prev_umid1) # arma_umid <- arma(umid,lag=list(ar=c(1,2,15)),include.intercept=T); arma_umid arima_umid <- arima(umid, order=c(15,0,0), fixed=c(NA,NA,0,0,0,0,0,0,0,0,0,0,0,0,NA,NA),include.mean=T) arima_umid prev_umid <- forecast(arima_umid, h = 12); prev_umid plot(prev_umid) # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAPÍTULO 10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # MODELOS SAZONAIS ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ################################################################################ # MODELOS SAZONAIS ESTOCÁSTICOS ################################################################################ ################################################################################ # Exemplo 10.8 Série IPI ################################################################################ IPI <- read.csv('C:\\R\\IPI.csv', head=T);IPI ipi <- IPI[1:180,2];ipi length(ipi) ts.plot(ipi, main="Precipitação de jan-66 a dez-97", xlab = "anos") summary(ipi) # ------------------------------------------------------------------------------ # Identificação do modelo preliminar # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(ipi, ci.type="ma") pacf(ipi) par(mfrow=c(2,2)) acf(ipi, lag.max = 36, main='d=0 e D=0', ci.type="ma") acf(diff(ipi), lag.max = 36, main='d=1 e D=0', ci.type="ma") acf(diff(ipi,lag=12), lag.max = 36, main='d=0 e D=1', ci.type="ma") acf(diff(diff(ipi),lag=12), lag.max = 36, main='d=1 e D=1', ci.type="ma") par(mfrow=c(1,2)) acf(diff(ipi,lag=12), lag.max = 36, main='d=0 e D=1', ci.type="ma") pacf(diff(ipi,lag=12), lag.max = 36, main='d=0 e D=1') acf(diff(diff(ipi),lag=12), lag.max = 36, main='d=1 e D=1', ci.type="ma") pacf(diff(diff(ipi),lag=12), lag.max = 36, main='d=1 e D=1') # Modelos preliminares: # m1: SARIMA(0,0,1)x(0,1,1)12 # m2: SARIMA(0,1,1)x(0,1,1)12 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar 10.2 (theta1, THETA1) # ------------------------------------------------------------------------------ m2 <- arima(ipi,order=c(0,1,1),seasonal = list(order = c(0, 1, 1), period = 12)) m2 # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar m2 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m2$resid, main='autocorrelações',ci.type="ma") pacf(m2$resid, main='autocorrelações parciais') # verifica-se valores altos de r4, phi4,4 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar m1 modificado incompleto (phi4, theta1, theta4 # THETA1 # ------------------------------------------------------------------------------ m22 <- arima(ipi,order=c(4,0,4),fixed=c(0,0,0,NA,NA,0,0,NA,NA),seasonal = list(order = c(0, 1, 1), period = 12), transform.pars=FALSE) m22 # testar a significância dos coeficientes estimados!!!!!!!!!! # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar modificado m22 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m22$resid, main='autocorrelações',ci.type="ma") pacf(m22$resid, main='autocorrelações parciais') # verifica-se valores altos de r2 e phi2,2, phi9,9 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar m111 modificado incompleto (phi2, phi4, phi9, # theta1, theta2, theta4, THETA1 # ------------------------------------------------------------------------------ m222 <- arima(ipi,order=c(9,0,4),fixed=c(0,NA,0,NA,0,0,0,0,NA,NA,NA,0,NA,NA), seasonal = list(order = c(0, 1, 1), period = 12), transform.pars=FALSE) m222 # testar a significância dos coeficientes estimados!!!!!!!!!! # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar modificado m222 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m222$resid, main='autocorrelações',ci.type="ma") pacf(m222$resid, main='autocorrelações parciais') # verifica-se valores altos de phi14,14 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar m111 modificado incompleto (phi2, phi4, phi9, # phi14, theta1, theta2, theta4, THETA1 # ------------------------------------------------------------------------------ m2222 <- arima(ipi,order=c(9,0,14),fixed=c(0,NA,0,NA,0,0,0,0,NA,NA,NA,0,NA,0,0,0,0,0,0,0,0,0,NA,NA),seasonal = list(order = c(0, 1, 1), period = 12),transform.pars=FALSE) m2222 # testar a significância dos coeficientes estimados!!!!!!!!!! # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar modificado m111 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m2222$resid, main='autocorrelações',ci.type="ma") pacf(m2222$resid, main='autocorrelações parciais') # O modelo está OK!!!!!!!!! # ------------------------------------------------------------------------------ # Previsão com base no modelo ajustado m111 # ------------------------------------------------------------------------------ require(forecast) prev_m2222 <- forecast(m2222, h = 6); prev_m2222 plot(prev_m2222) ############### & ################ # Modelos preliminares: # m1: SARIMA(0,0,1)x(0,1,1)12 # m2: SARIMA(0,1,1)x(0,1,1)12 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar 10.1 # ------------------------------------------------------------------------------ m1 <- arima(ipi,order=c(0,0,1),seasonal = list(order = c(0, 1, 1), period = 12)) m1 # Com constante # theta0 <- seq(1,length(ipi),1) # m1 <- arima(ipi,order=c(0,0,1),seasonal = list(order = c(0, 1, 1), period = 12), # xreg = theta0);m1 # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar m1 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m1$resid, main='autocorrelações',ci.type="ma") pacf(m1$resid, main='autocorrelações parciais') # verifica-se valores altos de r2, phi99 e phi14,14 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar m1 modificado incompleto # ------------------------------------------------------------------------------ m11 <- arima(ipi,order=c(14,0,2),fixed=c(0,0,0,0,0,0,0,0,NA,0,0,0,0,NA,NA,NA,NA),seasonal = list(order = c(0, 1, 1), period = 12), transform.pars=FALSE) m11 # Com cte # m11 <- arima(ipi,order=c(14,0,2),fixed=c(0,0,0,0,0,0,0,0,NA,0,0,0,0,NA,NA,NA,NA, # NA),seasonal = list(order = c(0, 1, 1), period = 12), xreg = theta0);m11 # testar a significância dos coeficientes estimados!!!!!!!!!! # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar modificado m11 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m11$resid, main='autocorrelações',ci.type="ma") pacf(m11$resid, main='autocorrelações parciais', ylab=' ') # verifica-se valores altos de r3 e r7 # ------------------------------------------------------------------------------ # Estimação do modelo preliminar m111 modificado incompleto # ------------------------------------------------------------------------------ m111 <- arima(ipi,order=c(14,0,7),fixed=c(0,0,0,0,0,0,0,0,NA,0,0,0,0,NA,NA,NA,NA,0,0,0,NA,NA),seasonal = list(order = c(0, 1, 1), period = 12),transform.pars=FALSE) m111 # m111 <- arima(ipi,order=c(14,0,2),fixed=c(0,0,NA,0,0,0,NA,0,NA,0,0,0,0,NA,NA, # NA,NA,NA),seasonal = list(order = c(0, 1, 1), period = 12),xreg = theta0) # m111 # testar a significância dos coeficientes estimados!!!!!!!!!! # ------------------------------------------------------------------------------ # Diagnóstico do ajuste do modelo preliminar modificado m111 # ------------------------------------------------------------------------------ par(mfrow=c(1,2)) acf(m111$resid, main='autocorrelações',ci.type="ma") pacf(m111$resid, main='autocorrelações parciais', ylab=' ') # O modelo está OK!!!!!!!!! # ------------------------------------------------------------------------------ # Previsão com base no modelo ajustado m111 # ------------------------------------------------------------------------------ require(forecast) prev_m111 <- forecast(m111, h = 6); prev_m111 plot(prev_m111)