| Rev. | 502367225f8f688864c2cb2f46280ee4a8167be4 |
|---|---|
| 大小 | 8,860 字节 |
| 时间 | 2016-03-12 00:24:41 |
| 作者 | Lorenzo Isella |
| Log Message | A code for nowcasting berd also using some interesting time series transformations (including the removal of NaN etc...). |
Forkrm(list=ls())
library(ggplot2)
library(grid)
library(reshape2)
library(RJSDMX)
require(gridExtra)
library(tikzDevice)
library(scales)
library(digest)
library(readr)
library(viridis)
library(tis)
library(forecast)
library(dplyr)
source("/home/lorenzo/myprojects-hg/R-codes/stat_lib.R")
year_cut <- 2000
options( scipen = 16 )
golden_ratio <- golden()
data1 <- getTimeSeries('EUROSTAT',
'gba_nabsfin07.A.NABS99.MIO_NAC.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
data_gba <- reshape_series(data1, year_cut)
data2 <- getTimeSeries('EUROSTAT',
'rd_e_gerdtot.A.TOTAL|BES.MIO_NAC.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
data_gerd <- reshape_series(data2, year_cut)
data_gba$indicator <- "GBAORD"
sel <- grep("TOTAL", data_gerd$indicator)
data_gerd$indicator[sel] <- "GERD"
data_gerd$indicator[-sel] <- "BERD"
data_tot <- rbind(data_gba, data_gerd)
data_tot$indicator <- as.factor(data_tot$indicator)
lbls <- levels(data_tot$indicator)
col_seq <- seq(length(lbls))
#col_seq[5] <- "brown"
title_exp <- paste("R&D Expenditure and Allocations", sep="")
gpl <- ggplot(data_tot, aes(x=year, y=value ,
colour=indicator,
shape=indicator
)) +
my_ggplot_theme(c(0.13, 0.8))+
theme(legend.position = 'right')+
geom_point(size=3) +
geom_line()+
facet_wrap( ~ country, nrow = 7 , scales = "free_y"
)+
scale_color_viridis("",labels=c("BERD", "GBAORD", "GERD"),
discrete=TRUE,na.value = "grey85"## , formatter=comma
)+
scale_shape_manual("", breaks=lbls,
labels=c("BERD", "GBAORD", "GERD"),
values=col_seq) +
## scale_colour_manual("Sectors of performance\n", breaks=lbls,
## labels=c("BERD", "GBAORD", "GERD"),
## values=col_seq) +
#scale_y_continuous(limits=c(0.3,0.9),breaks=seq(0.3, 0.9, by=0.3))+
scale_y_continuous(breaks=pretty_breaks(n=4))+
scale_x_continuous(breaks=pretty_breaks(n=4))+
labs(title=title_exp)+
theme(plot.title = element_text(lineheight=.8, size=24, face="bold", vjust=1))+
theme(legend.text = element_text(vjust=1,lineheight=1 ))+
theme(legend.title = element_text(vjust=1,lineheight=1, size=12 ))+
xlab("Year")+
ylab("Millions National Currency")
fname <- paste("gerd-berd-gbaord.pdf", sep="")
ggsave(fname, gpl, width=10*golden_ratio,height=10)
fname <- paste("_gerd-berd-gbaord.png", sep="")
ggsave(fname, gpl, width=10*golden_ratio,height=10)
#################################################
##Some examples about DE
data3 <- getTimeSeries('EUROSTAT',
'namq_10_a10_e.Q.THS_PER.M_N.NSA.EMP_DC.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
de_emp_mn <- data3$namq_10_a10_e.Q.THS_PER.M_N.NSA.EMP_DC.DE
data4 <- getTimeSeries('EUROSTAT',
'namq_10_a10_e.Q.THS_PER.C.NSA.EMP_DC.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
de_emp_c <- data4$namq_10_a10_e.Q.THS_PER.C.NSA.EMP_DC.DE
data5 <- getTimeSeries('EUROSTAT',
'namq_10_a10.Q.CP_MNAC.NSA.M_N.B1G.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
de_gva_mn <- data5$namq_10_a10.Q.CP_MNAC.NSA.M_N.B1G.DE
data6 <- getTimeSeries('EUROSTAT',
'namq_10_a10.Q.CP_MNAC.NSA.C.B1G.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
de_gva_c <- data6$namq_10_a10.Q.CP_MNAC.NSA.C.B1G.DE
data7 <- getTimeSeries('EUROSTAT',
'namq_10_gdp.Q.CP_MNAC.NSA.B1GQ|P3|P62|P61.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
data8 <- getTimeSeries('EUROSTAT',
'namq_10_a10.Q.CP_MNAC.NSA.B-E.B1G.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
data9 <- getTimeSeries('EUROSTAT',
'namq_10_a10_e.Q.THS_PER.B-E.NSA.EMP_DC.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
de_emp_ind <- data9$"namq_10_a10_e.Q.THS_PER.B-E.NSA.EMP_DC.DE"
de_emp_ind <- na.omit(as.ts(de_emp_ind))
de_emp_ind_agg <- aggregate(de_emp_ind, nfreq=1)
de_exp_ind <- data8$"namq_10_a10.Q.CP_MNAC.NSA.B-E.B1G.DE"
de_exp_ind <- na.omit(as.ts(de_exp_ind))
de_exp_ind_agg <- aggregate(de_exp_ind, nfreq=1)
de_gdp <- data7$namq_10_gdp.Q.CP_MNAC.NSA.B1GQ.DE
de_gdp <- na.omit(as.ts(de_gdp))
de_gdp_agg <- aggregate(de_gdp, nfreq=1)
de_exp_ser <- data7$namq_10_gdp.Q.CP_MNAC.NSA.P62.DE
de_exp_ser <- na.omit(as.ts(de_exp_ser))
de_exp_ser_agg <- aggregate(de_exp_ser, nfreq=1)
## de_exp_ser <- window(de_exp_ser, year_cut)
de_exp_goods <- data7$namq_10_gdp.Q.CP_MNAC.NSA.P61.DE
de_exp_goods <- na.omit(as.ts(de_exp_goods))
de_exp_goods_agg <- aggregate(de_exp_goods, nfreq=1)
de_gov_exp <- data7$namq_10_gdp.Q.CP_MNAC.NSA.P3.DE
de_gov_exp <- na.omit(as.ts(de_gov_exp))
de_gov_exp_agg <- aggregate(de_gov_exp, nfreq=1)
de_gbaord <- data1$gba_nabsfin07.A.NABS99.MIO_NAC.DE
de_gbaord <- na.omit(as.ts(de_gbaord))
de_gbaord <- window(de_gbaord, 1991)
de_berd <- data2$rd_e_gerdtot.A.BES.MIO_NAC.DE
de_berd <- na.omit(as.ts(de_berd))
de_berd <- window(de_berd, 1991)
de_emp_mn <- na.omit(as.ts(de_emp_mn))
de_emp_mn_agg <- aggregate(de_emp_mn, nfreq=1)
de_emp_c <- na.omit(as.ts(de_emp_c))
de_emp_c_agg <- aggregate(de_emp_c, nfreq=1)
de_gva_mn <- na.omit(as.ts(de_gva_mn))
de_gva_mn_agg <- aggregate(de_gva_mn, nfreq=1)
de_gva_c <- na.omit(as.ts(de_gva_c))
de_gva_c_agg <- aggregate(de_gva_c, nfreq=1)
predictors <- cbind(de_gdp_agg, de_gbaord , de_exp_ser_agg, de_exp_goods_agg,
de_emp_mn_agg, de_emp_c_agg, de_gva_mn_agg, de_gva_c_agg,
de_gov_exp_agg, de_emp_ind_agg, de_exp_ind_agg
)
## ## I can download the data directly as a data frame in a wide normalized format
## data7bis <- getTimeSeriesTable('EUROSTAT',
## 'namq_10_gdp.Q.CP_MNAC.NSA.B1GQ|P3|P62|P61.BE|BG|CZ|DK|DE|EE|IE|EL|ES|FR|IT|CY|LV|LT|LU|HU|MT|NL|AT|PL|PT|RO|SI|SK|FI|SE|UK|HR')
data_df <- as.data.frame(predictors)
data_df_red <- data_df[1:(nrow(data_df)-1), ]
data_df_red$berd <- de_berd
mylm <- lm(berd~ de_emp_mn_agg+ de_emp_c_agg +
de_gov_exp_agg ## +
## de_emp_ind_agg
, data=data_df_red)
berd_new <- predict(mylm, data_df)
df <- data.frame(berd = c(de_berd), year = c(time(de_berd)))
df2 <- data.frame(berd =c(berd_new), year=c(df$year, (max(df$year)+1)))
df$type <- "historical data"
df2$type <- "model"
df2$type[nrow(df2)] <- "nowcast"
df_tot <- rbind(df, df2)
df_tot$type <- as.factor(df_tot$type)
lbls <- levels(df_tot$type)
col_seq <- seq(length(lbls))
#col_seq[5] <- "brown"
title_exp <- paste("German BERD Nowcast", sep="")
gpl <- ggplot(df_tot, aes(x=year, y=berd ,
colour=type,
shape=type
)) +
my_ggplot_theme(c(0.13, 0.8))+
theme(legend.position = 'right')+
geom_point(size=3) +
geom_line()+
## facet_wrap( ~ country, nrow = 7 , scales = "free_y"
## )+
scale_color_viridis(NULL,labels=c("data", "model", "nowcast"),
discrete=TRUE,na.value = "grey85"## , formatter=comma
)+
scale_shape_manual(NULL, breaks=lbls,
labels=c("data", "model", "nowcast"),
values=col_seq) +
## scale_colour_manual("Sectors of performance\n", breaks=lbls,
## labels=c("BERD", "GBAORD", "GERD"),
## values=col_seq) +
#scale_y_continuous(limits=c(0.3,0.9),breaks=seq(0.3, 0.9, by=0.3))+
scale_y_continuous(breaks=pretty_breaks(n=4))+
scale_x_continuous(breaks=pretty_breaks(n=4))+
labs(title=title_exp)+
theme(plot.title = element_text(lineheight=.8, size=24, face="bold", vjust=1))+
theme(legend.text = element_text(vjust=1,lineheight=1 ))+
theme(legend.title = element_text(vjust=1,lineheight=1, size=12 ))+
xlab("Year")+
ylab("BERD (MIO NAC)")
fname <- paste("berd_nowcast.pdf", sep="")
ggsave(fname, gpl, width=5*golden_ratio,height=5)
fname <- paste("berd_nowcast.png", sep="")
ggsave(fname, gpl, width=5*golden_ratio,height=5)
###########################################################à
## other attempts do not yield much
data_df_red_log <- log(data_df_red)
data_df_log <- log(data_df)
mylm_log <- lm(berd~ de_emp_mn_agg+ de_emp_c_agg +de_gov_exp_agg , data=data_df_red_log)
berd_new_log <- exp(predict(mylm_log, data_df_log))
xreg <- select(data_df_red, de_emp_mn_agg, de_emp_c_agg, de_gov_exp_agg )
arima_mm <- auto.arima(data_df_red$berd,xreg=xreg, stepwise=FALSE, approximation=FALSE,max.order=9)
print("So far so good")