Uncategorized

Momento R do Dia: nascem mais Akiras do que Massakos no Japão?

nicesexratios
A gurizada nasce…mas tem sempre um pouco mais de demanda por soldadinhos relativamente à de bonecas?

Intuitivamente, qual deve ser a razão dos sexos em um país? Digo, você espera que o percentual de nascimentos entre homens e mulheres seja de 50% ou diferente de 50%?

Pois é. Agora, considere os dados do Japão (tabela 4 desta página). A série histórica começa em 1947 e vai até 2010. A medida: The sex ratio at birth is the ratio of males to 100 females. Eu, que não sou da área de saúde, fiquei meio confuso com a medida, mas, ao consultar uma fonte superficial, vi que o que eu imaginava faz sentido. Assim, dividi os dados por 100.

Será esta uma série estacionária? Veja um resumo da série.

sex_ratiojapan1

Pois é. Dá para ficar em dúvida, não dá? Vejamos o gráfico da série mais de perto. Sim, você já o viu lá no alto. Ambos foram gerados no R, mas não vou dar a dica do gráfico mais bonito (ahá!). Aquele lá no alto é facilmente construído no R, mas, desta vez, não vou te dar almoço tão barato (quase grátis). Das dicas anteriores, você poderá facilmente descobrir como fazer um gráfico apresentável em relatórios. Lembre-se: estudar faz bem e eu investi tempo aqui. Portanto, faça sua parte também. Mas vamos em frente. Olha o gráfico aí.

sexratio_japan2_correto

O exercício, então, é o de fazer o teste para procurar evidências de raiz unitária na série. Será que a razão dos sexos é um passeio aleatório? Usando o teste ADF com tendência (não parece o melhor teste, não é?), percebe-se que o coeficiente da tendência determinista não é significativo. Já o teste apenas com deslocamento (drift) retorna (saída parcial):

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)      0.3169 0.1274 2.487 0.0158 *
z.lag.1          -0.2995 0.1204 -2.487 0.0158 *
z.diff.lag       -0.5293 0.1071 -4.944 6.88e-06 ***

Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.004193 on 58 degrees of freedom
Multiple R-squared: 0.5307, Adjusted R-squared: 0.5145
F-statistic: 32.79 on 2 and 58 DF, p-value: 2.973e-10

Value of test-statistic is: -2.4866 3.1163

Critical values for test statistics:
1pct 5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1 6.70 4.71 3.86

Ou seja, a hipótese de raiz unitária não é rejeitada e nem a de que seja um passeio aleatório puro. Assim, talvez a variação da razão dos sexos seja estacionária. Fazer o teste não é difícil, mas é bom dar uma olhada no gráfico da mesma antes.

diff_sex

Ok, há uma variação marcada ali no meio da década dos 60, já vista no gráfico anterior. Basta olhar os dados para ver o que há. Usei o default do teste KPSS (o exercício, aqui, é fazer o aluno aprender os comandos, então, vamos variar) para o teste sem tendência determinista, mas com deslocamento.

O resultado?

#######################
# KPSS Unit Root Test #
#######################

Test is of type: mu with 3 lags.

Value of test-statistic is: 0.0539

Critical value for a significance level of:
10pct 5pct 2.5pct 1pct
critical values     0.347 0.463 0.574 0.739

Bem, o valor calculado é menor do que qualquer valor crítico e a hipótese nula do teste é a de estacionaridade. Bem, a partir daqui pode-se começar a brincadeira de se estimar o modelo ARIMA que explica a diff(jsex_ratio).

Claro, um exercício mais completo exige que o aluno faça ambos os testes para ambas as séries: em nível e em diferença. Obviamente, sem ler a parte teórica antes, a repetição pura e simples não o permitirá uma nota de aprovação na prova.

Mas vejamos o que nos resta.

diff_sex_arima

Os correlogramas não mentem: na função de autocorrelação acumulada, vê-se claramente que as duas primeiras autocorrelações são significativas (estatisticamente diferentes de zero, caso você não tenha entendido). No caso da função parcial, apenas a primeira autocorrelação se destaca. Agora, o que isso significa? Uma opção é que seja um ARIMA(1,1,0). Mas a série é curta e anual, o que pode nos esconder os padrões. Por que não um ARIMA (0,1,2)? Ou ainda um ARIMA (1,1,1)?

A metodologia Box-Jenkins implicaria testar estes modelos. Vou apresentar os resultados para um deles, o ARIMA (1,1,1). Os outros ficam por sua conta. A seguir, os resíduos, mostrando-nos um bom ajuste.

arima111Veja lá. Os correlogramas do resíduos não mostram nada fora dos intervalos de confiança (obviamente, a defasagem zero não conta (por que?)). Os p-valores para o teste de Ljung-Box nos mostram uma alta probabilidade que não-rejeição da hipótese nula do teste (qual é?).

O modelo:

> arima111<-Arima(jsex_ratio,order=c(1,1,1))
> summary(arima111)
Series: jsex_ratio
ARIMA(1,1,1)

Coefficients:
ar1          ma1
-0.5179   -0.3523
s.e. 0.1608    0.1886

sigma^2 estimated as 1.814e-05: log likelihood=250.07
AIC=-494.14 AICc=-493.72 BIC=-487.75

Training set error measures:
ME               RMSE            MAE                  MPE           MAPE      MASE
Training set 1.00896e-05 0.004227749 0.002987203 -0.0002905378 0.2817057 0.7698836

Pois é. Com estas medidas de acurácia de previsão, podemos comparar estes e outros modelos. Sobre a metodologia Box-Jenkins, bem, veremos depois. Para facilitar sua vida, eis o código usado aqui (com mais resultados do que os apresentados).

sex_ratio <- read.table(file = "clipboard", sep = "\t", header=TRUE)
head(sex_ratio)

japan_sex_ratio<-ts(sex_ratio$Sex.ratio.at.birth,start=c(1947),freq=1)

plot(japan_sex_ratio)

japan_sex_ratio_adjusted<-japan_sex_ratio/100
jsex_ratio<-japan_sex_ratio_adjusted
plot(jsex_ratio)

library(forecast)
tsdisplay(jsex_ratio)

library(urca)
summary(ur.df(jsex_ratio,type=c("trend")))
summary(ur.df(jsex_ratio,type=c("drift")))

plot(diff(jsex_ratio))

summary(ur.kpss(diff(jsex_ratio), type=c("mu")))

tsdisplay(diff(jsex_ratio))
arima111<-Arima(jsex_ratio,order=c(1,1,1))
summary(arima111)

# salvamos os resíduos e examinamos...
res<-residuals(arima111)
tsdisplay(res)

# fazemos o teste de Ljung-Box

Box.test(res, lag=10, fitdf=4, type="Ljung")

# ou poderíamos ver os residuos assim:
tsdiag(arima111)

# e a previsão...

plot(forecast(arima111,h=10)

Bem, pessoal, então, é isto. A razão dos sexos é um passeio aleatório puro e a parte estacionária pode ser modelada como um ARIMA (1,1,1). O exercício mostra o que mais você pode fazer. Primeiro, após ler sobre a metodologia, estimar os outros modelos e gerar resultados similares aos deste texto. Claro, em segundo lugar, compará-los. Em terceiro lugar, não se esquecer de olhar não apenas os gráficos, mas os testes de diagnóstico, bem como as medidas de acurácia de previsão.

Aqueles que curtem um teste sobre estas medidas de acurácia, busquem na documentação do R algo sobre o teste de Diebold-Mariano, por exemplo.

Concluindo, as mulheres devem continuar sendo alvo de disputas entre os garotos lá no Japão, não?

Um comentário em “Momento R do Dia: nascem mais Akiras do que Massakos no Japão?

Deixe um comentário

Preencha os seus dados abaixo ou clique em um ícone para log in:

Logotipo do WordPress.com

Você está comentando utilizando sua conta WordPress.com. Sair /  Alterar )

Foto do Google

Você está comentando utilizando sua conta Google. Sair /  Alterar )

Imagem do Twitter

Você está comentando utilizando sua conta Twitter. Sair /  Alterar )

Foto do Facebook

Você está comentando utilizando sua conta Facebook. Sair /  Alterar )

Conectando a %s