#Downloading needed libraries
library(tidyverse)
library(broom)
library(haven)
library(AER)
library(sandwich)
library(lmtest)
library(car)
library(stargazer)
library(dplyr)
library(ggplot2)
library(GGally)
library(psych)
library(corrplot)
library(readxl)
library(writexl)
#TS graphs
library(quantmod)
library(tseries)
library(forecast)
library(TSstudio)
library(reshape2)
library(plotly)
#Histograms
library(TSstudio)
library(reshape2)
library(tidyverse)
library(hrbrthemes)
library(viridis)
library(forcats)
#Downloading needed data
dat_1 = read_xls('C:/Users/Kolya/Documents/R/R files/Term_paper/TP_data/OCVED_FULL_05_10.xls')
length(unique(dat_1$...2)) #83
## [1] 85
dat_2 = read_xls('C:/Users/Kolya/Documents/R/R files/Term_paper/TP_data/OCVED_FULL_11_17.xls')
length(unique(dat_2$...2)) #87
## [1] 89
dat_extra = read_xlsx('C:/Users/Kolya/Documents/R/R files/Term_paper/TP_data/Classes.xlsx')
names(dat_1) <- c(c('industry', 'region', 'r.2005.1', 'r.2005.2', 'r.2005.3'), paste("r", sort(rep(c(2006:2010), 4)), rep(c(1,2,3,4),4), sep=".") )
dat_1 <- dat_1[4:length(dat_1$industry),]
names(dat_2) <- c(c('industry', 'region'), paste("r", sort(rep(c(2011:2016), 4)),rep(c(1,2,3,4),4), sep="."), 'y.2017')
dat_2 <- dat_2[4:length(dat_2$industry),]
dat <- merge(dat_1, dat_2, by = c('industry', 'region'), sort = FALSE, all.x = TRUE, all.y = TRUE)
unique(dat$region)
## [1] "Российская Федерация"
## [2] "Белгородская область"
## [3] "Брянская область"
## [4] "Владимирская область"
## [5] "Воронежская область"
## [6] "Ивановская область"
## [7] "Калужская область"
## [8] "Костромская область"
## [9] "Курская область"
## [10] "Липецкая область"
## [11] "Московская область"
## [12] "Орловская область"
## [13] "Рязанская область"
## [14] "Смоленская область"
## [15] "Тамбовская область"
## [16] "Тверская область"
## [17] "Тульская область"
## [18] "Ярославская область"
## [19] "Город Москва столица Российской Федерации город федерального значения"
## [20] "Республика Карелия"
## [21] "Республика Коми"
## [22] "Архангельская область"
## [23] "Ненецкий автономный округ (Архангельская область)"
## [24] "Вологодская область"
## [25] "Калининградская область"
## [26] "Ленинградская область"
## [27] "Мурманская область"
## [28] "Новгородская область"
## [29] "Псковская область"
## [30] "Город Санкт-Петербург город федерального значения"
## [31] "Республика Адыгея (Адыгея)"
## [32] "Республика Калмыкия"
## [33] "Краснодарский край"
## [34] "Астраханская область"
## [35] "Волгоградская область"
## [36] "Ростовская область"
## [37] "Республика Дагестан"
## [38] "Республика Ингушетия"
## [39] "Кабардино-Балкарская Республика"
## [40] "Карачаево-Черкесская Республика"
## [41] "Республика Северная Осетия-Алания"
## [42] "Чеченская Республика"
## [43] "Ставропольский край"
## [44] "Республика Башкортостан"
## [45] "Республика Марий Эл"
## [46] "Республика Мордовия"
## [47] "Республика Татарстан (Татарстан)"
## [48] "Удмуртская Республика"
## [49] "Чувашская Республика - Чувашия"
## [50] "Пермский край"
## [51] "Кировская область"
## [52] "Нижегородская область"
## [53] "Оренбургская область"
## [54] "Пензенская область"
## [55] "Самарская область"
## [56] "Саратовская область"
## [57] "Ульяновская область"
## [58] "Курганская область"
## [59] "Свердловская область"
## [60] "Тюменская область"
## [61] "Ханты-Мансийский автономный округ - Югра (Тюменская область)"
## [62] "Ямало-Ненецкий автономный округ (Тюменская область)"
## [63] "Челябинская область"
## [64] "Республика Алтай"
## [65] "Республика Тыва"
## [66] "Республика Хакасия"
## [67] "Алтайский край"
## [68] "Красноярский край"
## [69] "Иркутская область"
## [70] "Кемеровская область - Кузбасс"
## [71] "Новосибирская область"
## [72] "Омская область"
## [73] "Томская область"
## [74] "Республика Бурятия"
## [75] "Забайкальский край"
## [76] "Республика Саха (Якутия)"
## [77] "Камчатский край"
## [78] "Приморский край"
## [79] "Хабаровский край"
## [80] "Амурская область"
## [81] "Магаданская область"
## [82] "Сахалинская область"
## [83] "Еврейская автономная область"
## [84] "Чукотский автономный округ"
## [85] "Архангельская область (без АО)"
## [86] "Республика Крым"
## [87] "Город федерального значения Севастополь"
## [88] "Тюменская область (без АО)"
col_for_sum <- c()
tmp <- dat[, c(3:50)] %>% mutate_if(is.character,as.numeric)
for(year in c(2005:2016)){
year <- toString(year)
subyear <- tmp[, grepl(year, names(tmp))]
tmp[paste('r_m',year, sep = '_')] <- rowMeans(subyear)
col_for_sum <- append(col_for_sum, paste('r_m',year, sep = '_'))
}
tmp$r_m_2017 <- tmp$y.2017
tmp <- tmp[, c(49:61)]
tmp <- cbind(dat[,c(1,2)], tmp)
tmp <- tmp[c(tmp$region != 'Российская Федерация' &
tmp$region != 'Ненецкий автономный округ (Архангельская область)' &
tmp$region != 'Ямало-Ненецкий автономный округ (Тюменская область)' &
tmp$region != 'Архангельская область (без АО)' &
tmp$region != 'Тюменская область (без АО)' &
tmp$region != 'Республика Крым' &
tmp$region != 'Город федерального значения Севастополь' &
tmp$region != 'Чеченская Республика' &
tmp$region != 'Ханты-Мансийский автономный округ - Югра (Тюменская область)'),]
#Replacing names with shorter ones
tmp$region <-
replace(tmp$region,
tmp$region == "Город Москва столица Российской Федерации город федерального значения",
'Москва')
tmp$region <-
replace(tmp$region, tmp$region == "Город Санкт-Петербург город федерального значения",
"Санкт-Петербург")
unique(tmp$region)
## [1] "Белгородская область" "Брянская область"
## [3] "Владимирская область" "Воронежская область"
## [5] "Ивановская область" "Калужская область"
## [7] "Костромская область" "Курская область"
## [9] "Липецкая область" "Московская область"
## [11] "Орловская область" "Рязанская область"
## [13] "Смоленская область" "Тамбовская область"
## [15] "Тверская область" "Тульская область"
## [17] "Ярославская область" "Москва"
## [19] "Республика Карелия" "Республика Коми"
## [21] "Архангельская область" "Вологодская область"
## [23] "Калининградская область" "Ленинградская область"
## [25] "Мурманская область" "Новгородская область"
## [27] "Псковская область" "Санкт-Петербург"
## [29] "Республика Адыгея (Адыгея)" "Республика Калмыкия"
## [31] "Краснодарский край" "Астраханская область"
## [33] "Волгоградская область" "Ростовская область"
## [35] "Республика Дагестан" "Республика Ингушетия"
## [37] "Кабардино-Балкарская Республика" "Карачаево-Черкесская Республика"
## [39] "Республика Северная Осетия-Алания" "Ставропольский край"
## [41] "Республика Башкортостан" "Республика Марий Эл"
## [43] "Республика Мордовия" "Республика Татарстан (Татарстан)"
## [45] "Удмуртская Республика" "Чувашская Республика - Чувашия"
## [47] "Пермский край" "Кировская область"
## [49] "Нижегородская область" "Оренбургская область"
## [51] "Пензенская область" "Самарская область"
## [53] "Саратовская область" "Ульяновская область"
## [55] "Курганская область" "Свердловская область"
## [57] "Тюменская область" "Челябинская область"
## [59] "Республика Алтай" "Республика Тыва"
## [61] "Республика Хакасия" "Алтайский край"
## [63] "Красноярский край" "Иркутская область"
## [65] "Кемеровская область - Кузбасс" "Новосибирская область"
## [67] "Омская область" "Томская область"
## [69] "Республика Бурятия" "Забайкальский край"
## [71] "Республика Саха (Якутия)" "Камчатский край"
## [73] "Приморский край" "Хабаровский край"
## [75] "Амурская область" "Магаданская область"
## [77] "Сахалинская область" "Еврейская автономная область"
## [79] "Чукотский автономный округ"
dat_extra <- dat_extra[,c(2,5,7:9)]
tmp$industry <- gsub("(РАЗДЕЛ\\s+.\\s+|Подраздел\\s+..\\s+)", "", tmp$industry)
tmp$industry <- toupper(tmp$industry)
tmp$industry <- str_squish(tmp$industry)
dat_extra$Name <- toupper(dat_extra$Name)
dat_extra$Name <- str_squish(dat_extra$Name)
tmp_1 <- merge(tmp, dat_extra, by.x = 'industry', by.y = 'Name', sort = TRUE)
tmp_1[,3:15][is.na(tmp_1[, 3:15])]<- 0.0000001
#Loosing is caused by bad classification in our data
length(setdiff(unique(tmp$industry), unique(dat_extra$Name)))
## [1] 294
#It should be enough for the folowing analysis
length(unique(tmp_1$industry))
## [1] 1453
#needed for calculation of shares pi
col_for_sum <- c(col_for_sum, 'r_m_2017')
# level 1
set_level_1 <- subset(tmp_1, is.na(tmp_1$SubKod1) == F & is.na(tmp_1$SubKod2) == T & is.na(tmp_1$SubKod3) == T)
sum_by_region <- aggregate(set_level_1[,col_for_sum], by=list(set_level_1$region), sum, na.rm=T)
names(sum_by_region)[1] <- 'region'
set_level_1 <- merge(set_level_1, sum_by_region, by=c('region'))
set_level_5 <- subset(tmp_1, is.na(tmp_1$SubKod3) == F)
set_level_5 <- merge(set_level_5, sum_by_region, by=c('region'))
Pi <- set_level_5[, 3:15]/set_level_5[,20:32]
Pi$region <- set_level_5$region
Pi$industry <- set_level_5$industry
Pi$SubKod1 <- set_level_5$SubKod1
Pi$SubKod2 <- set_level_5$SubKod2
Pi$SubKod3 <- set_level_5$SubKod3
names(Pi)[1:13] <- paste('Pi', c(2005:2017), sep = '_')
#Check
Pi[Pi$region == 'Алтайский край' & Pi$industry == 'РАСПРЕДЕЛЕНИЕ ВОДЫ', ]
#set_for_Pg <- merge(Pi, tmp_1[, c(1,2,17:19)], by=c('region', 'industry'))
Pg <- aggregate(Pi[,1:13], by=list(Pi$region, Pi$SubKod1, Pi$SubKod2), sum, na.rm=T)
names(Pg) <- c('region', 'SubKod1', 'SubKod2', paste('Pg', c(2005:2017), sep = '_'))
#Check
Pg[Pg$SubKod2 == '00' & Pg$region=='Алтайский край', ]
set_for_Hg <- merge(Pg, Pi[,c(1:14,16,17, 18)], by=c('region', 'SubKod1', 'SubKod2'))
set_for_Hg <- cbind(set_for_Hg, set_for_Hg[, 17:29]/set_for_Hg[, 4:16]*log(1/(set_for_Hg[, 17:29]/set_for_Hg[,4:16]), base = 2))
names(set_for_Hg)[31:43] <- paste('Hg', c(2005:2017), sep = '_')
Hg <- aggregate(set_for_Hg[,31:43], by=list(set_for_Hg$region, set_for_Hg$SubKod1, set_for_Hg$SubKod2), sum, na.rm=T)
names(Hg)[1:3] <- c('region', 'SubKod1', 'SubKod2')
RV
set_for_Rv <- merge(Pg, Hg, by=c('region', 'SubKod1', 'SubKod2'))
set_for_Rv <- cbind(set_for_Rv, set_for_Rv[,4:16]*set_for_Rv[,17:29])
names(set_for_Rv)[30:42] <- paste('Rv_part', c(2005:2017), sep = '_')
RV <- aggregate(set_for_Rv[,30:42], by=list(set_for_Rv$region), sum, na.rm=T)
names(RV)[1] <- c('region')
names(RV)[2:14] <- paste('Rv', c(2005:2017), sep = '_')
aggregate(RV$Rv_2005, by=list(RV$region), sum)
RV
UV
Pg_for_UV <- Pg
Pg_for_UV <- cbind(Pg_for_UV, Pg_for_UV[,4:16]*log(1/Pg_for_UV[,4:16],base = 2))
names(Pg_for_UV)[17:29] <- paste('Pg_log', c(2005:2017), sep = '_')
UV <- aggregate(Pg_for_UV[,17:29], by=list(Pg_for_UV$region), sum, na.rm=T)
names(UV)[1:14] <- c('region', paste('Uv', c(2005:2017), sep = '_'))
UV
V
Set_for_Var <- Pi
Set_for_Var <- cbind(Set_for_Var, Pi[,1:13]*log(1/Pi[,1:13], base = 2))
Variety <- aggregate(Set_for_Var[,19:31], by=list(Set_for_Var$region), sum, na.rm=T)
names(Variety)[1:14] <- c('region', paste('V', c(2005:2017), sep = '_'))
Variety
Descriptive statistics
RV_stat <- RV %>% describe()
RV_stat <- RV_stat[2:14, ]
RV_stat <- colMeans(RV_stat[,c(3:4, 8:9)])
RV_stat
## mean sd min max
## 0.57412252 0.39532358 0.05489331 2.94027954
UV_stat <- UV %>% describe()
UV_stat <- UV_stat[2:14, ]
UV_stat <- colMeans(UV_stat[,c(3:4, 8:9)])
UV_stat
## mean sd min max
## 2.7446030 0.6276677 1.2006333 4.4163654
V_stat <- Variety %>% describe()
V_stat <- V_stat[2:14, ]
V_stat <- colMeans(V_stat[,c(3:4, 8:9)])
V_stat
## mean sd min max
## 3.318726 0.892582 1.285390 6.416477
Statistics <- t(data.frame(RV_stat, UV_stat, V_stat))
row.names(Statistics) <- c('Related Variety', 'Unrelated Variety', 'Variety')
stargazer(Statistics, title="Descriptive statistics", type = 'text', align=T, digits=3)
##
## Descriptive statistics
## =========================================
## mean sd min max
## -----------------------------------------
## Related Variety 0.574 0.395 0.055 2.940
## Unrelated Variety 2.745 0.628 1.201 4.416
## Variety 3.319 0.893 1.285 6.416
## -----------------------------------------
Time Series analysis
# RV TS
RV_series <- ts(colMeans(RV[2:14]), start = 2005, end = 2017)
# chartSeries(RV_series)
ts_plot(RV_series,
title = "mean Related Variety",
Xtitle = "Time",
Ytitle = "",
color = "green",
Xgrid = TRUE,
Ygrid = TRUE) %>%
layout(paper_bgcolor = "black",
plot_bgcolor = "black",
font = list(color = "white"),
yaxis = list(linecolor = "#6b6b6b",
zerolinecolor = "#6b6b6b",
gridcolor= "#444444"),
xaxis = list(linecolor = "#6b6b6b",
zerolinecolor = "#6b6b6b",
gridcolor= "#444444"))
# UV TS
UV_series <- ts(colMeans(UV[2:14]), start = 2005, end = 2017)
#chartSeries(UV_series)
ts_plot(UV_series,
title = "mean Unrelated Variety",
Xtitle = "Time",
Ytitle = "",
color = "green",
Xgrid = TRUE,
Ygrid = TRUE) %>%
layout(paper_bgcolor = "black",
plot_bgcolor = "black",
font = list(color = "white"),
yaxis = list(linecolor = "#6b6b6b",
zerolinecolor = "#6b6b6b",
gridcolor= "#444444"),
xaxis = list(linecolor = "#6b6b6b",
zerolinecolor = "#6b6b6b",
gridcolor= "#444444"))
# Variety TS
V_series <- ts(colMeans(Variety[2:14]), start = 2005, end = 2017)
#chartSeries(V_series)
ts_plot(V_series,
title = "mean Variety",
Xtitle = "Time",
Ytitle = "",
color = "green",
Xgrid = TRUE,
Ygrid = TRUE) %>%
layout(paper_bgcolor = "black",
plot_bgcolor = "black",
font = list(color = "white"),
yaxis = list(linecolor = "#6b6b6b",
zerolinecolor = "#6b6b6b",
gridcolor= "#444444"),
xaxis = list(linecolor = "#6b6b6b",
zerolinecolor = "#6b6b6b",
gridcolor= "#444444"))
meltdf <- melt(Variety,id="region")
meltdf$region <- as.factor(meltdf$region)
regions <- c(Variety$region[1:10])
df <- meltdf[meltdf$region %in% regions, ]
ggplot(df, aes(x=variable, y=value, colour=region, group=region)) +
geom_line() +
geom_point() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.1, hjust=2))
#Histogram for all years
ggplot(data = meltdf, aes(x=value, fill=variable)) +
geom_histogram(alpha=0.4, position = 'identity', color="#e9ecef", binwidth = 0.4) +
scale_fill_viridis(discrete=TRUE) +
scale_color_viridis(discrete=TRUE) +
theme_ipsum()
# ggsave('my_plot_2.png')
#Histograms for every year
ggplot(data = meltdf, aes(x=value, fill=variable)) +
geom_histogram(alpha=0.4, position = 'identity', color="#e9ecef", binwidth = 0.4) +
scale_fill_viridis(discrete=TRUE) +
scale_color_viridis(discrete=TRUE) +
theme_ipsum() +
theme(
legend.position="none",
panel.spacing = unit(0.1, "lines"),
strip.text.x = element_text(size = 8),
strip.text.y = element_text(size = 8)
) +
facet_wrap(~variable)
write_xlsx(Variety, 'C:/Users/Kolya/Documents/New_studies/Term_paper/Data/Raw data/V_index.xlsx')
write_xlsx(RV, 'C:/Users/Kolya/Documents/New_studies/Term_paper/Data/Raw data/RV_index.xlsx')
write_xlsx(UV, 'C:/Users/Kolya/Documents/New_studies/Term_paper/Data/Raw data/UV_index.xlsx')