Ciclo de Negócios - PIB SP

Para a execução do código são necessários os pacotes instalados com o seguinte comando: install.packages(c("mFilter", "BCDating", "ggplot2", "ggfortify", "gridExtra", "zoo"))

1. Carregamento de pacotes necessários

library(mFilter)
library(BCDating)
library(ggplot2)
library(ggfortify)
library(gridExtra)
library(zoo)

O pacote zoo é necessário para a conversão de datas de objetos ts.

2. Carregamento os dados

Dados do PIB mensal paulista da Fundação SEADE.

O código abaixo carrega os dados do PIB mensal no arquivo PIBSPDESMENSAL.txt e transforma num objeto ts do R.

PIBSPDESMENSAL <- read.table("PIBSPDESMENSAL.txt", header = 1)
pibspdes.ts <- ts(PIBSPDESMENSAL[,1], start = c(2002, 1), freq = 12)
pibspdes.ts

A Time Series: 18 × 12

	Jan	Feb	Mar	Apr	May	Jun	Jul	Aug	Sep	Oct	Nov	Dec
2002	73.33478	73.06779	73.64706	71.76099	70.93057	70.97745	69.59280	69.97219	71.06672	71.05543	68.92039	70.46595
2003	70.25221	69.91517	70.06311	70.19920	70.42013	70.38047	71.24865	71.89990	73.40967	72.77346	72.00712	72.95622
2004	72.99799	74.67012	74.62521	74.61258	75.67317	75.66943	75.87418	76.60303	75.89974	76.03392	77.05789	76.92885
2005	77.58925	77.85808	77.65135	78.17482	78.46795	78.62313	78.94936	78.86570	78.64828	80.32922	80.18617	80.88479
2006	80.74956	80.28435	80.90182	81.09761	80.63397	81.48458	81.74595	81.56003	82.69094	83.25254	83.89853	84.22039
2007	84.80598	83.96821	86.27099	87.28801	87.49391	87.84886	89.03589	89.68530	90.02649	89.50613	90.84226	90.72636
2008	88.66001	92.33906	92.11330	93.19064	94.07928	95.24395	94.17783	96.54803	94.11640	93.53354	88.83879	88.33278
2009	90.14950	88.96092	91.34150	91.45888	92.15164	93.11597	93.56399	94.96909	96.08541	97.84143	96.93722	99.22027
2010	98.66773	99.52808	99.44595	98.74035	100.17127	98.90729	101.86711	100.11137	100.89978	101.18079	101.65303	102.72070
2011	103.69311	103.49000	101.83190	104.99832	103.64968	105.39131	103.80004	104.68971	103.27766	104.76684	104.62375	103.97065
2012	102.62789	104.61862	103.81797	104.43381	104.47445	106.18359	107.49942	105.41123	106.18121	105.63915	106.69299	106.12409
2013	106.50347	107.00528	108.48911	108.01314	108.95287	108.50516	110.26984	110.09879	109.66959	109.34650	106.88631	108.79727
2014	108.66700	111.18375	107.20996	107.18646	105.36554	105.26512	105.98165	105.31337	106.57447	106.50425	105.27362	106.09808
2015	105.90576	104.82888	103.57051	103.12076	102.77488	101.02893	100.96674	100.63923	101.15105	99.52379	100.43213	99.93101
2016	98.32905	100.09144	100.13205	99.67750	99.66865	99.79342	98.27034	98.67453	98.07267	98.28473	97.74856	99.34015
2017	100.39278	97.88011	99.43144	100.00270	101.19828	101.58511	101.61454	102.10672	101.66481	102.44880	102.51445	101.71462
2018	101.36487	101.89256	102.37728	99.18121	103.60518	102.47171	102.84375	102.33969	101.81237	102.74005	103.65017	102.84418
2019	103.84362	102.08969	103.76201

3. Aplicação da função de datação de ciclos de negócios

Para obter a datação de ciclo de negócios é aplicada a função BBQ. A função BBQ, do pacote BCDating recebe o objeto de série temporal com valores transformados em logaritmo e multiplicado por 100, retornando um objeto com as datas de início e fim dos ciclos. O objeto retornado é atribuído à variável dc.

Para plotar os dados com o ggplot2 é conveniente que estejam no formato de dataframe. Por isso constrói-se um dataframe atribuindo-o à variável business.cycle.

dc <- BBQ(
    log(pibspdes.ts)*100,
    name="Dating Business Cycles of LOG PIBSP DESSAZONALIZADO")

a <- c(as.Date(dc@y)[1], as.Date(dc@y)[dc@peaks])
b <- as.Date(dc@y)[dc@troughs]
business.cycle <- data.frame(Topo=a, Vale=b)
business.cycle

A data.frame: 16 × 2

Topo	Vale
<date>	<date>
2002-01-01	2002-11-01
2003-09-01	2003-11-01
2004-02-01	2004-04-01
2005-07-01	2005-09-01
2005-12-01	2006-02-01
2007-11-01	2008-01-01
2008-08-01	2008-12-01
2010-02-01	2010-04-01
2011-01-01	2011-03-01
2011-06-01	2012-01-01
2013-07-01	2013-11-01
2014-02-01	2014-06-01
2014-09-01	2014-11-01
2015-09-01	2016-01-01
2016-03-01	2016-11-01
2017-11-01	2018-04-01

4. Aplicação o Filtro HP

Para obter os componentes Trend e Cycle do PIB paulista utiliza-se a função hpfilter, que recebe como argumento a série temporal do PIB.

Semelhante ao caso anterior, é conveniente que os dados estejam na forma de dataframe, nesse caso atribuído à variável pibspdes.ts.

pibsphp <- hpfilter(pibspdes.ts, freq = 14400, type = c("lambda"), drift = FALSE)
a <- cbind(pibsphp$x, pibsphp$trend, pibsphp$cycle)
pibspdes.ts <- data.frame(as.matrix(a), as.Date(as.yearmon(time(a))))
colnames(pibspdes.ts) <- c("PIB", "Trend", "Cycle", "Date")
pibspdes.ts

A data.frame: 207 × 4

PIB	Trend	Cycle	Date
<dbl>	<dbl>	<dbl>	<date>
73.33478	70.05551	3.27926554	2002-01-01
73.06779	70.15481	2.91297373	2002-02-01
73.64706	70.25435	3.39270943	2002-03-01
71.76099	70.35454	1.40645848	2002-04-01
70.93057	70.45605	0.47451734	2002-05-01
70.97745	70.55965	0.41780225	2002-06-01
69.59280	70.66613	-1.07333531	2002-07-01
69.97219	70.77632	-0.80413185	2002-08-01
71.06672	70.89097	0.17574440	2002-09-01
71.05543	71.01078	0.04465433	2002-10-01
68.92039	71.13644	-2.21605346	2002-11-01
70.46595	71.26868	-0.80272306	2002-12-01
70.25221	71.40804	-1.15582981	2003-01-01
69.91517	71.55503	-1.63986164	2003-02-01
70.06311	71.71007	-1.64696014	2003-03-01
70.19920	71.87346	-1.67426184	2003-04-01
70.42013	72.04540	-1.62527056	2003-05-01
70.38047	72.22596	-1.84549141	2003-06-01
71.24865	72.41511	-1.16645522	2003-07-01
71.89990	72.61268	-0.71277378	2003-08-01
73.40967	72.81842	0.59125408	2003-09-01
72.77346	73.03204	-0.25857609	2003-10-01
72.00712	73.25328	-1.24616043	2003-11-01
72.95622	73.48187	-0.52565286	2003-12-01
72.99799	73.71746	-0.71946359	2004-01-01
74.67012	73.95964	0.71048286	2004-02-01
74.62521	74.20797	0.41724415	2004-03-01
74.61258	74.46205	0.15052586	2004-04-01
75.67317	74.72152	0.95164762	2004-05-01
75.66943	74.98602	0.68341286	2004-06-01
...	...	...	...
98.28473	100.4987	-2.21394779	2016-10-01
97.74856	100.4694	-2.72082994	2016-11-01
99.34015	100.4547	-1.11458449	2016-12-01
100.39278	100.4541	-0.06127485	2017-01-01
97.88011	100.4666	-2.58650647	2017-02-01
99.43144	100.4917	-1.06021854	2017-03-01
100.00270	100.5283	-0.52558145	2017-04-01
101.19828	100.5755	0.62281080	2017-05-01
101.58511	100.6322	0.95290819	2017-06-01
101.61454	100.6975	0.91706532	2017-07-01
102.10672	100.7704	1.33633192	2017-08-01
101.66481	100.8501	0.81475222	2017-09-01
102.44880	100.9357	1.51307866	2017-10-01
102.51445	101.0267	1.48777030	2017-11-01
101.71462	101.1223	0.59230373	2017-12-01
101.36487	101.2221	0.14273022	2018-01-01
101.89256	101.3257	0.56688669	2018-02-01
102.37728	101.4325	0.94479678	2018-03-01
99.18121	101.5421	-2.36093715	2018-04-01
103.60518	101.6543	1.95086139	2018-05-01
102.47171	101.7685	0.70322306	2018-06-01
102.84375	101.8843	0.95946876	2018-07-01
102.33969	102.0014	0.33832128	2018-08-01
101.81237	102.1195	-0.30713322	2018-09-01
102.74005	102.2384	0.50160562	2018-10-01
103.65017	102.3579	1.29224780	2018-11-01
102.84418	102.4777	0.36644926	2018-12-01
103.84362	102.5977	1.24589173	2019-01-01
102.08969	102.7178	-0.62812602	2019-02-01
103.76201	102.8380	0.92405246	2019-03-01

5. Plotando os dados

Para plotar os dados utiliza-se o pacote ggplot2. Abaixo são salvas algumas variáveis que serão usadas nas configurações dos gráficos.

plot.limits <- as.Date(c("2010-12-31", "2019-03-02"))
plot.date_breaks <- "3 month"
plot.date_minor_breaks <- "1 month"
xscale.date.size <- 6

Para plotar o gráfico usa-se a função ggplot com a série temporal pibspdes.ts e a coluna Date como o eixo x.

Para que os dados apareçam no gráfico adiciona-se (+) à saída da função ggplot os chamados 'geomas', dizendo como devem ser apresentados os dados.

g1 <- ggplot(pibspdes.ts, aes(x = Date)) +
    geom_line(aes(y = PIB), color = "black") +
    geom_line(aes(y = Trend), color = "red") + 
    ylim(95, 115) +
    geom_rect(
        data = business.cycle,
        inherit.aes = FALSE,
        aes(xmin = Topo, xmax = Vale, ymin = -Inf, ymax = Inf),
        alpha = 0.2) +
    theme_bw() +
    scale_x_date(
        date_labels = "%Y-%m",
        date_breaks = plot.date_breaks,
        date_minor_breaks = plot.date_minor_breaks,
        limits = plot.limits) +
    ylab("PIB SP dessazonalizado") +
    ggtitle("PIB SP e Ciclo de Negocios") +
    theme(
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        plot.title = element_text(lineheight = .8, face = "bold"))
g1

Warning message:
"Removed 108 rows containing missing values (geom_path)."Warning message:
"Removed 108 rows containing missing values (geom_path)."Warning message:
"Removed 8 rows containing missing values (geom_rect)."

g2 <- ggplot(pibspdes.ts, aes(x = Date, y = Cycle)) +
    geom_col() +
    geom_rect(
        data = business.cycle,
        inherit.aes=FALSE,
        aes(xmin = Topo, xmax = Vale, ymin = -Inf, ymax = Inf),
        alpha = 0.2) +
    theme_bw() +
    scale_x_date(
        date_labels = "%Y-%m",
        date_breaks = plot.date_breaks,
        date_minor_breaks = plot.date_minor_breaks,
        limits = plot.limits) +
    ylab("Filtro HP - Ciclo") +
    theme(
        axis.text.x = element_text(
            angle = 90, 
            vjust = 0.25, 
            size = xscale.date.size))
g2

Warning message:
"Removed 108 rows containing missing values (position_stack)."Warning message:
"Removed 2 rows containing missing values (geom_col)."Warning message:
"Removed 8 rows containing missing values (geom_rect)."

Para que os gráficos fiquem numa única figura usa-se a função grid.arrange.

g <- grid.arrange(g1, g2, ncol = 1, heights = c(2, 1))
g

Warning message:
"Removed 108 rows containing missing values (geom_path)."Warning message:
"Removed 108 rows containing missing values (geom_path)."Warning message:
"Removed 8 rows containing missing values (geom_rect)."Warning message:
"Removed 108 rows containing missing values (position_stack)."Warning message:
"Removed 2 rows containing missing values (geom_col)."Warning message:
"Removed 8 rows containing missing values (geom_rect)."


TableGrob (2 x 1) "arrange": 2 grobs
  z     cells    name           grob
1 1 (1-1,1-1) arrange gtable[layout]
2 2 (2-2,1-1) arrange gtable[layout]

ggsave(file = "plot.png", plot = g, width = 18, height = 9, dpi = 300)