# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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=100, list(ar=c(0.3))) par(mfrow=c(1,3)) acf.ar <- acf(ar.sim) acf.ar <- acf(ar.sim, ci.type="ma") pacf(ar.sim) plot(acf.ar,ci.type="ma",ci=0.90) plot(acf.ar,ci.type="ma",ci=0.95) plot(acf.ar,ci.type="ma",ci=0.99) ################################################################################ # Exemplo 6.3 - série ICV ################################################################################ # ICV <- read.csv("D:\\Users\\User\\Documents\\icv.txt", head=T) ICV <- read.csv("G:\\Meu Drive\\esalq\\pos_graduacao\\diciplina_series_temporais\\TS\\2022\\icv.txt", head=T) ICV icv <- ICV[,1] icv plot(icv, main="ICV de jan-70 a jun-79", lwd = "2", col = "51", xlab = "anos", ylim = c(10,850)) 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") # ------------------------------------------------------------------------------ # Quantas diferen?as? # ------------------------------------------------------------------------------ # veja a FAC - 1o lag: ~ -0.5 # veja a FACP - mudan?a de sinais dos valores da FACP # e quanto a vari?ncia? # Vari?ncia da 1a diferen?a vard1logicv <- var(d1logicv) vard1logicv # Vari?ncia da 2a diferen?a vard2logicv <- var(d2logicv) vard2logicv # ------------------------------------------------------------------------------ # Teste de Raiz Unit?ria (Dickey-Fuller) # H0: h? raiz unit?ria # ------------------------------------------------------------------------------ require(urca) ?ur.df logicv <- log(icv) logicv1 <- ts(logicv) # teste com intercepto, sem tend?ncia linear e sem dummies sazonais teste1 <- ur.df(logicv1, type=c("drift"), lags=0) teste2 <- ur.df(d1logicv1, type=c("drift"), lags=0) summary(teste1) summary(teste2) # H0 not rejected # teste com intercepto, sem tend?ncia linear e sem dummies sazonais d1logicv <- diff(logicv) d1logicv1 <- ts(d1logicv) teste2 <- ADF.test (wts=d1logicv1,selectlags=list(Pmax=1),itsd=c(1,0,c(0)), regvar=0) teste2 # H0 rejected at 1% significancy, thus d=1 # ------------------------------------------------------------------------------ # Testar se a m?dia ? estatisticamente diferente de zero (t-ratio) # ------------------------------------------------------------------------------ var_d1logicv <- var(d1logicv) r1 <- acf(d1logicv) r1$acf[2] varW_bar <- var_d1logicv*(1+r1$acf[2])/(113*(1-r1$acf[2])) varW_bar t_ratio <- mean(d1logicv)/(varW_bar)^0.5;t_ratio # H0 ? rejeitado # ------------------------------------------------------------------------------ # 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.csv2("G:\\Meu Drive\\esalq\\pos_graduacao\\diciplina_series_temporais\\TS\\2022\\atmosfera.csv", head=T, sep = ",", dec = ".") UMID names(UMID) umid <- UMID [1:358,3];umid 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') # funcao BoxCox.ar e BoxCox help.search("boxcox") require(forecast) require(TSA) ?BoxCox ?BoxCox.ar umid1 <- as.ts(umid) umid1 plot(umid1) BoxCox.ar(umid1) par(mfrow=c(2,1)) plot(BoxCox(umid1,2)) plot(umid1) # ------------------------------------------------------------------------------ # 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 # ------------------------------------------------------------------------------ # Teste de Raiz Unit?ria (Dickey-Fuller) # H0: h? raiz unit?ria # ------------------------------------------------------------------------------ require(uroot) ADF.test umid1 <- ts(umid) # teste com intercepto, sem tend?ncia linear e sem dummies sazonais teste <- ADF.test (wts=umid1,selectlags=list(Pmax=1),itsd=c(0,0,c(0)), regvar=0) teste # H0 rejected # ------------------------------------------------------------------------------ # Testar se a m?dia ? estatisticamente diferente de zero (t-ratio) # H0: E(W) = 0. # ------------------------------------------------------------------------------ var_umid <- var(umid) r1 <- acf(umid) r1$acf[2] varW_bar <- var_umid*(1+r1$acf[2])/(357*(1-r1$acf[2])) varW_bar t_ratio <- mean(umid)/(varW_bar)^0.5;t_ratio # H0 ? rejeitado # ------------------------------------------------------------------------------ # 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.csv2("C:\\Users\\vitor\\OneDrive\\Documentos\\poluicao.csv", head=T, dec=".") #CO CO <- read.csv2("G:\\Meu Drive\\esalq\\pos_graduacao\\diciplina_series_temporais\\TS\\2022\\poluicao.csv", head=T, sep = ",", dec = ".") 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(co) # ------------------------------------------------------------------------------ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CAP?TULO 7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ################################################################################ # ESTIMA??O DE MODELOS ARIMA(p,d,q) ################################################################################ # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ################################################################################ # Exemplo 7.1 EMV-Procedimento Condicional ################################################################################ #Zt <- c(150,147,143,148,153,149,155,162,170,172) #Zt # Ajuste de um modelo IMA(1,1) #ima <- arima(Zt,order=c(0,1,1));ima # Obs: o m?nimo de S*(theta) ocorre ao redor de theta = 0.4. # Vari?ncia do par?metro #var_theta <- (1-0.4^2)/9 #ima$var.coef ################################################################################ # 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 # ------------------------------------------------------------------------------ # 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 # ------------------------------------------------------------------------------ # 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) theta0 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). # ARIMAX 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!! # ------------------------------------------------------------------------------ # C?lculo da constante # ------------------------------------------------------------------------------ theta0 <- 0.0237*(1-0.5024); theta0 # ------------------------------------------------------------------------------ # C?lculo do erro padr?o # ------------------------------------------------------------------------------ ep <- (arima_icv1$var.coef)^0.5;ep # ------------------------------------------------------------------------------ # Signific?ncia dos par?metros # ------------------------------------------------------------------------------ # phi1 0.5073/0.0822 # intercept 0.0212/0.0018 # ------------------------------------------------------------------------------ # Outra forma de estimar modelos ARIMA(p,d,q) # ------------------------------------------------------------------------------ #require(tseries) #arma_icv <- arma(d1logicv,order=c(1,0), include.intercept=T); arma_icv # nessa fun??o o resultado "intercept" ? de fato a constante! ################################################################################ # 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 # ------------------------------------------------------------------------------ # 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 # ------------------------------------------------------------------------------ # 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