Uncategorized

Momento R do Dia – Penn World Tables 8.0 e a Propensão Média a Consumir de cada país

Fazia tempo, heim? Mas hoje um aluno me perguntou sobre algo que falei no início do semestre (é…alguém não prestou atenção…): consumo. Mais especificamente, sobre a propensão média a consumir.

Eu acredito que o bom aluno não tem apenas que perguntar, mas também pesquisar e pensar nas implicações do que pergunta. Para ajudar, vou fazer um cálculo simples e rápido que ilustra dois pacotes do R, um que apenas tem a base de dados (“pwt8”) e o outro que aplica uma função para calcular a propensão média a consumir média para cada país.

Explico: a base de dados tem 127 países e cada qual tem variáveis para o período 1950-2011. Nem todas as variáveis estão, entretanto, registradas, para todos os países. Assim, por exemplo, a propensão média a consumir que calculo para o Brasil a seguir é a média das propensões médias deste período (uma propensão média “média” a consumir, por assim dizer). Já para Angola, que não tem os dados para todo o período, eu também calculo esta mesma variável, mas excluindo os períodos em que não há valores para uma, outra ou ambas as variáveis.

O código abaixo possui as explicações mais ou menos detalhadas.

# como sempre...você só instala uma vez. Ou seja, depois não precisa mais destas duas linhas de
# comando. O pacote plyr tem o famoso conjunto de comandos "apply" (lapply, sapply, etc).

install.packages("pwt8")
install.packages("plyr")

# carregamos o pacote, carregamos a base de dados do pacote.
library(pwt8)
data("pwt8.0")
names(pwt8.0)
# como exemplo, eis o pib per capita calculado para todos os anos, todos os países
pwt8.0$cap <- pwt8.0$rgdpe / pwt8.0$pop

# agora calculo a propensao media a conumir. Primeiro, calculamos o consumo, pela multiplicacao
# da participação do consumo sobre o PIB vezes o PIB. Veja a documentação de pwt8 para ver
# a definição de csh_c e cgdpe. Note que, na primeira linha, gero o consumo e, em seguida,
# na segunda linha, a Propensão Média a Consumir (C/Y) para cada ano. 

pwt8.0$cons<-(pwt8.0$csh_c)*pwt8.0$cgdpe
pwt8.0$pmec<-pwt8.0$cons / pwt8.0$cgdpe

# agora vamos criar um data.frame com a base e vamos chamá-la de "dta.
# Para ter uma idéia, usamos "head" para ver as primeiras linhas do arquivo.

dta <- data.frame(pwt8.0)
head(dta)

# tudo parece ok. Repare que, como indexei "cons" e "pmec" pela pwt8.0 lá em cima, agora elas
# fazem parte da base de dados pwt8.0 que, nestes comandos acima, transformei em uma estrutura
# de data.frame (não precisava, mas é útil fazê-lo para outras coisas...).
# Agora carregamos os comandos de plyr

library(plyr)

# vamos usar "ddply" (leia a documentação do pacote ou veja as dicas mais adiante no post)
# repare que o comando utiliza a base "dta", organiza por código de país ("isocode"),
# pede um sumário, por país da variável "pmec_media" que é a media das pmec calculadas
# anteriormente, ignorando os anos em que houve algum dado faltante ("na.rm=TRUE"). 

ddply(dta, ~ isocode, summarize, pmec_media = mean(pmec, na.rm=TRUE))

Agora, vejamos o resultado. Para ver mais detalhes, veja a documentação de pwt8 e plyr. Outras dicas são esta e esta.

isocode pmec_media
1 AGO 0.2698527
2 ALB 0.5749997
3 ARG 0.7098281
4 ARM 0.7371243
5 ATG 0.5383594
6 AUS 0.5989834
7 AUT 0.6242973
8 AZE 0.4728226
9 BDI 0.8484496
10 BEL 0.5939900
11 BEN 0.6626455
12 BFA 0.7358716
13 BGD 0.8368961
14 BGR 0.5891358
15 BHR 0.5313741
16 BHS 0.6305308
17 BIH 0.8930009
18 BLR 0.5099687
19 BLZ 0.7110097
20 BMU 0.9430771
21 BOL 0.6384951
22 BRA 0.6557681
23 BRB 0.9176295
24 BRN 0.1783046
25 BTN 0.4300709
26 BWA 0.5208991
27 CAF 0.6800754
28 CAN 0.6121375
29 CHE 0.5868197
30 CHL 0.6855206
31 CHN 0.5494199
32 CIV 0.7488742
33 CMR 0.7527335
34 COD 0.6749305
35 COG 0.4382671
36 COL 0.7529843
37 COM 0.5984154
38 CPV 0.6978079
39 CRI 0.8358405
40 CYP 0.5574062
41 CZE 0.5002645
42 DEU 0.5809886
43 DJI 0.6269397
44 DMA 0.6015103
45 DNK 0.5624745
46 DOM 0.7318125
47 ECU 0.6132193
48 EGY 0.6695550
49 ESP 0.6522298
50 EST 0.5297069
51 ETH 0.7811232
52 FIN 0.5024147
53 FJI 0.6368589
54 FRA 0.5974393
55 GAB 0.4527526
56 GBR 0.6328932
57 GEO 0.7129032
58 GHA 0.7278015
59 GIN 0.7263120
60 GMB 0.8160562
61 GNB 0.8856748
62 GNQ 0.4649407
63 GRC 0.6865936
64 GRD 0.7845137
65 GTM 0.8610200
66 HKG 0.5691085
67 HND 0.7746579
68 HRV 0.5953630
69 HUN 0.5360931
70 IDN 0.5866937
71 IND 0.7005359
72 IRL 0.6487952
73 IRN 0.4866217
74 IRQ 0.3241603
75 ISL 0.5204941
76 ISR 0.5020997
77 ITA 0.5891751
78 JAM 0.7392437
79 JOR 0.7085971
80 JPN 0.5370400
81 KAZ 0.5876354
82 KEN 0.8001836
83 KGZ 0.6266033
84 KHM 0.8267202
85 KNA 0.6442140
86 KOR 0.5968393
87 KWT 0.2970611
88 LAO 0.5390379
89 LBN 0.8995980
90 LBR 0.8636384
91 LCA 0.8534209
92 LKA 0.6632201
93 LSO 1.2219911
94 LTU 0.6099957
95 LUX 0.5703499
96 LVA 0.5631248
97 MAC 0.3930121
98 MAR 0.6934141
99 MDA 0.7019827
100 MDG 0.8074133
101 MDV 0.1832236
102 MEX 0.7255086
103 MKD 0.7062246
104 MLI 0.7376391
105 MLT 1.3556602
106 MNE 0.5621934
107 MNG 0.5627013
108 MOZ 0.9719093
109 MRT 0.6465644
110 MUS 0.7040063
111 MWI 0.7647157
112 MYS 0.5356852
113 NAM 0.6724800
114 NER 0.5999834
115 NGA 0.6010883
116 NLD 0.5312290
117 NOR 0.4914107
118 NPL 0.7126890
119 NZL 0.6175504
120 OMN 0.3364292
121 PAK 0.7486567
122 PAN 0.6806284
123 PER 0.7052688
124 PHL 0.6864944
125 POL 0.5937320
126 PRT 0.6926333
127 PRY 0.6934128
128 QAT 0.2613939
129 ROU 0.5638809
130 RUS 0.4876037
131 RWA 0.7233076
132 SAU 0.3245511
133 SDN 0.9571723
134 SEN 0.7633639
135 SGP 0.4474048
136 SLE 0.8428216
137 SLV 1.2248291
138 SRB 0.6809770
139 STP 1.1072268
140 SUR 0.6335978
141 SVK 0.5602101
142 SVN 0.5883688
143 SWE 0.5505727
144 SWZ 0.7613371
145 SYR 0.6650626
146 TCD 0.2787636
147 TGO 0.7384517
148 THA 0.5818592
149 TJK 0.5880398
150 TKM 0.3464478
151 TTO 0.6275485
152 TUN 0.6485752
153 TUR 0.6332986
154 TWN 0.5346673
155 TZA 0.4372538
156 UGA 0.6687535
157 UKR 0.5639214
158 URY 0.7252260
159 USA 0.6825835
160 UZB 0.3869788
161 VCT 0.7430303
162 VEN 0.4057366
163 VNM 0.7126157
164 YEM 0.6690849
165 ZAF 0.6422902
166 ZMB 0.7610378
167 ZWE 0.6292824

Bonito, heim? Repare que um ou dois países deveriam ter seus dados analisados com cautela e não é difícil entender o porquê.

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?