Some Convergence and Mobility Facts

library(tidyverse)
library(ggrepel)
library(extrafont)
library(viridis)
library(summarytools)
library(xtable)
library(knitr)
options(prompt="R> ", digits=2)

1 Sector ISIC18LC-stch

1.1 Sigma and Beta Convergence Analysis

# Data import
stch18 <- read_csv(" ISIC18LC-stch/inputs/dataMobility.csv")
stch18 <- stch18 %>% 
  mutate(id = as.factor(id))

# Compute average annual growth rates (use approximations)
stch18 <- stch18 %>% 
  mutate(growth1 = ((ln_y - ln_x)/6)*100,
         growth2 =  ((y/x)^(1/6)-1)*100)
stch18

# Descriptive statistics of continous variables
stch18Growth <- stch18 %>% 
  select(x, y, ln_x, ln_y, growth1, growth2, r_x, r_y, ln_r_x, ln_r_y) 
kable(descr(stch18Growth))

	x	y	ln_x	ln_y	growth1	growth2	r_x	r_y	ln_r_x	ln_r_y
Mean	3.01	3.19	0.76	0.79	0.44	2.07	1.00	1.00	-0.34	-0.37
Std.Dev	2.98	4.26	0.80	0.73	17.86	19.11	0.99	1.33	0.80	0.73
Min	0.35	0.63	-1.04	-0.46	-38.47	-31.93	0.12	0.20	-2.14	-1.63
Median	2.19	1.87	0.78	0.62	0.59	0.60	0.73	0.58	-0.32	-0.54
Max	15.69	22.76	2.75	3.13	49.63	64.26	5.22	7.12	1.65	1.96
MAD	1.40	0.78	0.70	0.46	12.10	11.69	0.47	0.24	0.70	0.46
IQR	1.68	1.16	0.85	0.61	17.28	16.86	0.56	0.36	0.85	0.61
CV	1.01	0.75	0.96	1.09	0.02	0.11	1.01	0.75	-0.42	-0.51
Skewness	2.29	3.42	0.29	1.44	0.37	1.05	2.29	3.42	0.29	1.44
SE.Skewness	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33
Kurtosis	5.54	11.80	0.03	2.11	0.64	1.67	5.54	11.80	0.03	2.11
N.Valid	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00
Pct.Valid	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

# Create long panel and recode year
stch18Panel <- stch18 %>%
  gather(year, ln_tfp, ln_x:ln_y)

# Recode year variable
stch18Panel$year <- fct_recode(stch18Panel$year,
                               "2003" = "ln_x",
                               "2009" = "ln_y")

# Plot basic Boxplot
stch18Panel %>% 
  ggplot(aes(x = year, y = ln_tfp, fill = year)) +
  geom_boxplot(outlier.colour="black", 
               outlier.shape=16,
               outlier.size=1,
               notch=TRUE) +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position = 'none') +
  theme(text=element_text(size=10, family="Palatino"))

# Beta convergence scatterplot 
stch18Growth %>% 
  ggplot(aes(x = ln_x, y = growth1)) +
  geom_point(size = 4) + 
  geom_smooth(size = 2) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of TFP in 2003",
       y = "TFP Growth")

1.2 Mobility Analysis

# Select and rename variables
stch18 <- stch18 %>% 
  select(id, ln_r_x, ln_r_y) %>% 
  rename(initial = ln_r_x, final = ln_r_y)

# Identify forward and backward mobility
stch18 <- stch18 %>% 
  mutate(mobility0 = as.factor(final > initial))

# Identify mobility across quadrants
stch18 <- stch18 %>% 
  mutate(mobility = as.factor(case_when(final > 0 & initial > 0 ~ 1, 
                              final > 0 & initial < 0 ~ 2,
                              final < 0 & initial < 0 ~ 3,
                              final < 0 & initial > 0 ~ 4,
                              TRUE ~ 0))
                              )
# Identify observation in the first quadrant (uppper right corner)
stch18 <- stch18 %>% 
  mutate(mobility1 = as.factor(final > 0 & initial > 0))
  
# Identify observation in the second quadrant (lower right corner)
stch18 <- stch18 %>% 
  mutate(mobility2 = as.factor(final > 0 & initial < 0))
# Identify observation in the third quadrant (lower left corner)
stch18 <- stch18 %>% 
  mutate(mobility3 = as.factor(final < 0 & initial < 0))
# Identify observation in the third quadrant (upper left corner)
stch18 <- stch18 %>% 
  mutate(mobility4 = as.factor(final < 0 & initial > 0))

# Print table
stch18

# Create scatterplot 
stch18 %>% 
ggplot(aes(x = final , y = initial, color = mobility0, shape = mobility)) +
  geom_point(size = 4) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Create scatterplot with bivarate density
stch18 %>% 
ggplot(aes(x = final , y = initial)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_gradientn(colors = c("#FFEDA0", "#FEB24C", "#F03B20"))+
  geom_point(size = 2) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Descriptive statistics for continuous variables
kable(descr(stch18))

	initial	final
Mean	-0.34	-0.37
Std.Dev	0.80	0.73
Min	-2.14	-1.63
Median	-0.32	-0.54
Max	1.65	1.96
MAD	0.70	0.46
IQR	0.85	0.61
CV	-0.42	-0.51
Skewness	0.29	1.44
SE.Skewness	0.33	0.33
Kurtosis	0.03	2.11
N.Valid	53.00	53.00
Pct.Valid	100.00	100.00

# Frequency table for forward and backward mobility
kable(freq(stch18$mobility0))

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
FALSE	27	51	51	51	51
TRUE	26	49	100	49	100
	0	NA	NA	0	100
Total	53	100	100	100	100

# Frequency table for mobility across quadrants
freq(stch18$mobility)

Frequencies   
mobility     
Data frame: stch18   
Type: Factor (unordered)   

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          1      4      7.55           7.55      7.55           7.55
          2      8     15.09          22.64     15.09          22.64
          3     32     60.38          83.02     60.38          83.02
          4      9     16.98         100.00     16.98         100.00
       <NA>      0                               0.00         100.00
      Total     53    100.00         100.00    100.00         100.00

# Simple mobility matrix
kable(with(stch18, ctable(mobility0, mobility, prop = "r")) )

	1	2	3	4	Total
FALSE	4	0	14	9	27
TRUE	0	8	18	0	26
Total	4	8	32	9	53

	1	2	3	4	Total
FALSE	0.15	0.00	0.52	0.33	1
TRUE	0.00	0.31	0.69	0.00	1
Total	0.08	0.15	0.60	0.17	1

2 Sector ISIC18LC-lp

2.1 Sigma and Beta Convergence Analysis

# Data import
lp18 <- read_csv(" ISIC18LC-lp/inputs/dataMobility.csv")
lp18 <- lp18 %>% 
  mutate(id = as.factor(id))

# Compute average annual growth rates (use approximations)
lp18 <- lp18 %>% 
  mutate(growth1 = ((ln_y - ln_x)/6)*100,
         growth2 =  ((y/x)^(1/6)-1)*100)
lp18

# Descriptive statistics of continous variables
lp18Growth <- lp18 %>% 
  select(x, y, ln_x, ln_y, growth1, growth2, r_x, r_y, ln_r_x, ln_r_y) 
kable(descr(lp18Growth))

	x	y	ln_x	ln_y	growth1	growth2	r_x	r_y	ln_r_x	ln_r_y
Mean	3.20	3.15	0.81	0.79	-0.31	1.29	1.00	1.00	-0.35	-0.36
Std.Dev	3.24	3.78	0.81	0.75	17.85	18.66	1.01	1.20	0.81	0.75
Min	0.44	0.53	-0.81	-0.63	-43.22	-35.09	0.14	0.17	-1.98	-1.78
Median	2.18	1.95	0.78	0.67	0.83	0.84	0.68	0.62	-0.38	-0.48
Max	16.86	20.82	2.82	3.04	47.85	61.36	5.27	6.61	1.66	1.89
MAD	1.44	0.91	0.78	0.54	12.32	12.95	0.45	0.29	0.78	0.54
IQR	2.16	1.13	1.02	0.61	22.23	21.65	0.68	0.36	1.02	0.61
CV	0.99	0.83	1.00	1.06	-0.02	0.07	0.99	0.83	-0.43	-0.48
Skewness	2.23	3.01	0.43	1.10	0.22	0.86	2.23	3.01	0.43	1.10
SE.Skewness	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33
Kurtosis	5.18	9.43	-0.23	1.05	0.35	1.27	5.18	9.43	-0.23	1.05
N.Valid	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00
Pct.Valid	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

# Create long panel and recode year
lp18Panel <- lp18 %>%
  gather(year, ln_tfp, ln_x:ln_y)

# Recode year variable
lp18Panel$year <- fct_recode(lp18Panel$year,
                               "2003" = "ln_x",
                               "2009" = "ln_y")

# Plot basic Boxplot
lp18Panel %>% 
  ggplot(aes(x = year, y = ln_tfp, fill = year)) +
  geom_boxplot(outlier.colour="black", 
               outlier.shape=16,
               outlier.size=1,
               notch=TRUE) +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position = 'none') +
  theme(text=element_text(size=10, family="Palatino"))

# Beta convergence scatterplot 
lp18Growth %>% 
  ggplot(aes(x = ln_x, y = growth1)) +
  geom_point(size = 4) + 
  geom_smooth(size = 2) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of TFP in 2003",
       y = "TFP Growth")

2.2 Mobility Analysis

# Select and rename variables
lp18 <- lp18 %>% 
  select(id, ln_r_x, ln_r_y) %>% 
  rename(initial = ln_r_x, final = ln_r_y)

# Identify forward and backward mobility
lp18 <- lp18 %>% 
  mutate(mobility0 = as.factor(final > initial))

# Identify mobility across quadrants
lp18 <- lp18 %>% 
  mutate(mobility = as.factor(case_when(final > 0 & initial > 0 ~ 1, 
                              final > 0 & initial < 0 ~ 2,
                              final < 0 & initial < 0 ~ 3,
                              final < 0 & initial > 0 ~ 4,
                              TRUE ~ 0))
                              )
# Identify observation in the first quadrant (uppper right corner)
lp18 <- lp18 %>% 
  mutate(mobility1 = as.factor(final > 0 & initial > 0))
  
# Identify observation in the second quadrant (lower right corner)
lp18 <- lp18 %>% 
  mutate(mobility2 = as.factor(final > 0 & initial < 0))
# Identify observation in the third quadrant (lower left corner)
lp18 <- lp18 %>% 
  mutate(mobility3 = as.factor(final < 0 & initial < 0))
# Identify observation in the third quadrant (upper left corner)
lp18 <- lp18 %>% 
  mutate(mobility4 = as.factor(final < 0 & initial > 0))

# Print table
lp18

# Create scatterplot 
lp18 %>% 
ggplot(aes(x = final , y = initial, color = mobility0, shape = mobility)) +
  geom_point(size = 4) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Create scatterplot with bivarate density
lp18 %>% 
ggplot(aes(x = final , y = initial)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_gradientn(colors = c("#FFEDA0", "#FEB24C", "#F03B20"))+
  geom_point(size = 2) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Descriptive statistics for continuous variables
kable(descr(lp18))

	initial	final
Mean	-0.35	-0.36
Std.Dev	0.81	0.75
Min	-1.98	-1.78
Median	-0.38	-0.48
Max	1.66	1.89
MAD	0.78	0.54
IQR	1.02	0.61
CV	-0.43	-0.48
Skewness	0.43	1.10
SE.Skewness	0.33	0.33
Kurtosis	-0.23	1.05
N.Valid	53.00	53.00
Pct.Valid	100.00	100.00

# Frequency table for forward and backward mobility
kable(freq(lp18$mobility0))

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
FALSE	24	45	45	45	45
TRUE	29	55	100	55	100
	0	NA	NA	0	100
Total	53	100	100	100	100

# Frequency table for mobility across quadrants
freq(lp18$mobility)

Frequencies   
mobility     
Data frame: lp18   
Type: Factor (unordered)   

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          1      4      7.55           7.55      7.55           7.55
          2      7     13.21          20.75     13.21          20.75
          3     31     58.49          79.25     58.49          79.25
          4     11     20.75         100.00     20.75         100.00
       <NA>      0                               0.00         100.00
      Total     53    100.00         100.00    100.00         100.00

# Simple mobility matrix
kable(with(lp18, ctable(mobility0, mobility, prop = "r")) )

	1	2	3	4	Total
FALSE	2	0	11	11	24
TRUE	2	7	20	0	29
Total	4	7	31	11	53

	1	2	3	4	Total
FALSE	0.08	0.00	0.46	0.46	1
TRUE	0.07	0.24	0.69	0.00	1
Total	0.08	0.13	0.58	0.21	1

3 Sector ISIC18LC-acf

3.1 Sigma and Beta Convergence Analysis

# Data import
acf18 <- read_csv(" ISIC18LC-acf/inputs/dataMobility.csv")
acf18 <- acf18 %>% 
  mutate(id = as.factor(id))

# Compute average annual growth rates (use approximations)
acf18 <- acf18 %>% 
  mutate(growth1 = ((ln_y - ln_x)/6)*100,
         growth2 =  ((y/x)^(1/6)-1)*100)
acf18

# Descriptive statistics of continous variables
acf18Growth <- acf18 %>% 
  select(x, y, ln_x, ln_y, growth1, growth2, r_x, r_y, ln_r_x, ln_r_y) 
kable(descr(acf18Growth))

	x	y	ln_x	ln_y	growth1	growth2	r_x	r_y	ln_r_x	ln_r_y
Mean	0.57	0.74	-0.94	-0.81	2.04	4.36	1.00	1.00	-0.37	-0.51
Std.Dev	0.57	1.25	0.87	0.85	21.11	23.13	1.00	1.68	0.87	0.85
Min	0.03	0.09	-3.41	-2.45	-54.56	-42.05	0.06	0.12	-2.85	-2.15
Median	0.40	0.35	-0.92	-1.04	2.14	2.16	0.70	0.48	-0.36	-0.74
Max	2.60	7.11	0.95	1.96	53.41	70.58	4.55	9.59	1.52	2.26
MAD	0.28	0.19	0.72	0.65	14.67	14.89	0.50	0.25	0.72	0.65
IQR	0.38	0.35	0.98	0.83	20.79	20.52	0.67	0.47	0.98	0.83
CV	1.00	0.59	-1.07	-0.96	0.10	0.19	1.00	0.59	-0.43	-0.60
Skewness	2.02	3.92	-0.03	1.26	0.28	1.04	2.02	3.92	-0.03	1.26
SE.Skewness	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33	0.33
Kurtosis	3.52	15.57	0.31	1.95	0.67	1.22	3.52	15.57	0.31	1.95
N.Valid	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00	53.00
Pct.Valid	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

# Create long panel and recode year
acf18Panel <- acf18 %>%
  gather(year, ln_tfp, ln_x:ln_y)

# Recode year variable
acf18Panel$year <- fct_recode(acf18Panel$year,
                               "2003" = "ln_x",
                               "2009" = "ln_y")

# Plot basic Boxplot
acf18Panel %>% 
  ggplot(aes(x = year, y = ln_tfp, fill = year)) +
  geom_boxplot(outlier.colour="black", 
               outlier.shape=16,
               outlier.size=1,
               notch=TRUE) +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position = 'none') +
  theme(text=element_text(size=10, family="Palatino"))

# Beta convergence scatterplot 
acf18Growth %>% 
  ggplot(aes(x = ln_x, y = growth1)) +
  geom_point(size = 4) + 
  geom_smooth(size = 2) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of TFP in 2003",
       y = "TFP Growth")

3.2 Mobility Analysis

# Select and rename variables
acf18 <- acf18 %>% 
  select(id, ln_r_x, ln_r_y) %>% 
  rename(initial = ln_r_x, final = ln_r_y)

# Identify forward and backward mobility
acf18 <- acf18 %>% 
  mutate(mobility0 = as.factor(final > initial))

# Identify mobility across quadrants
acf18 <- acf18 %>% 
  mutate(mobility = as.factor(case_when(final > 0 & initial > 0 ~ 1, 
                              final > 0 & initial < 0 ~ 2,
                              final < 0 & initial < 0 ~ 3,
                              final < 0 & initial > 0 ~ 4,
                              TRUE ~ 0))
                              )
# Identify observation in the first quadrant (uppper right corner)
acf18 <- acf18 %>% 
  mutate(mobility1 = as.factor(final > 0 & initial > 0))
  
# Identify observation in the second quadrant (lower right corner)
acf18 <- acf18 %>% 
  mutate(mobility2 = as.factor(final > 0 & initial < 0))
# Identify observation in the third quadrant (lower left corner)
acf18 <- acf18 %>% 
  mutate(mobility3 = as.factor(final < 0 & initial < 0))
# Identify observation in the third quadrant (upper left corner)
acf18 <- acf18 %>% 
  mutate(mobility4 = as.factor(final < 0 & initial > 0))

# Print table
acf18

# Create scatterplot 
acf18 %>% 
ggplot(aes(x = final , y = initial, color = mobility0, shape = mobility)) +
  geom_point(size = 4) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Create scatterplot with bivarate density
acf18 %>% 
ggplot(aes(x = final , y = initial)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_gradientn(colors = c("#FFEDA0", "#FEB24C", "#F03B20"))+
  geom_point(size = 2) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Descriptive statistics for continuous variables
kable(descr(acf18))

	initial	final
Mean	-0.37	-0.51
Std.Dev	0.87	0.85
Min	-2.85	-2.15
Median	-0.36	-0.74
Max	1.52	2.26
MAD	0.72	0.65
IQR	0.98	0.83
CV	-0.43	-0.60
Skewness	-0.03	1.26
SE.Skewness	0.33	0.33
Kurtosis	0.31	1.95
N.Valid	53.00	53.00
Pct.Valid	100.00	100.00

# Frequency table for forward and backward mobility
kable(freq(acf18$mobility0))

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
FALSE	31	58	58	58	58
TRUE	22	42	100	42	100
	0	NA	NA	0	100
Total	53	100	100	100	100

# Frequency table for mobility across quadrants
freq(acf18$mobility)

Frequencies   
mobility     
Data frame: acf18   
Type: Factor (unordered)   

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          1      3      5.66           5.66      5.66           5.66
          2      8     15.09          20.75     15.09          20.75
          3     29     54.72          75.47     54.72          75.47
          4     13     24.53         100.00     24.53         100.00
       <NA>      0                               0.00         100.00
      Total     53    100.00         100.00    100.00         100.00

# Simple mobility matrix
kable(with(acf18, ctable(mobility0, mobility, prop = "r")) )

	1	2	3	4	Total
FALSE	1	0	17	13	31
TRUE	2	8	12	0	22
Total	3	8	29	13	53

	1	2	3	4	Total
FALSE	0.03	0.00	0.55	0.42	1
TRUE	0.09	0.36	0.55	0.00	1
Total	0.06	0.15	0.55	0.25	1

4 Sector ISIC36LC-stch

4.1 Sigma and Beta Convergence Analysis

# Data import
stch36 <- read_csv(" ISIC36LC-stch/inputs/dataMobility.csv")
stch36 <- stch36 %>% 
  mutate(id = as.factor(id))

# Compute average annual growth rates (use approximations)
stch36 <- stch36 %>% 
  mutate(growth1 = ((ln_y - ln_x)/6)*100,
         growth2 =  ((y/x)^(1/6)-1)*100)
stch36

# Descriptive statistics of continous variables
stch36Growth <- stch36 %>% 
  select(x, y, ln_x, ln_y, growth1, growth2, r_x, r_y, ln_r_x, ln_r_y) 
kable(descr(stch36Growth))

	x	y	ln_x	ln_y	growth1	growth2	r_x	r_y	ln_r_x	ln_r_y
Mean	3.0	3.66	0.88	0.98	1.66	2.85	1.00	1.00	-0.21	-0.32
Std.Dev	2.2	4.03	0.65	0.73	14.95	16.58	0.73	1.10	0.65	0.73
Min	0.5	0.67	-0.69	-0.40	-24.90	-22.04	0.17	0.18	-1.78	-1.70
Median	2.4	2.43	0.87	0.89	-0.44	-0.44	0.80	0.66	-0.22	-0.41
Max	12.7	25.98	2.54	3.26	51.25	66.95	4.26	7.09	1.45	1.96
MAD	1.6	1.44	0.72	0.67	12.05	11.58	0.55	0.39	0.72	0.67
IQR	2.3	2.27	0.93	0.90	15.52	15.39	0.78	0.62	0.93	0.90
CV	1.4	0.91	1.35	1.34	0.11	0.17	1.37	0.91	-0.32	-0.44
Skewness	2.1	3.35	0.04	0.82	0.89	1.42	2.09	3.35	0.04	0.82
SE.Skewness	0.3	0.30	0.30	0.30	0.30	0.30	0.30	0.30	0.30	0.30
Kurtosis	5.7	13.72	-0.14	0.56	0.95	2.68	5.71	13.72	-0.14	0.56
N.Valid	62.0	62.00	62.00	62.00	62.00	62.00	62.00	62.00	62.00	62.00
Pct.Valid	100.0	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

# Create long panel and recode year
stch36Panel <- stch36 %>%
  gather(year, ln_tfp, ln_x:ln_y)

# Recode year variable
stch36Panel$year <- fct_recode(stch36Panel$year,
                               "2003" = "ln_x",
                               "2009" = "ln_y")

# Plot basic Boxplot
stch36Panel %>% 
  ggplot(aes(x = year, y = ln_tfp, fill = year)) +
  geom_boxplot(outlier.colour="black", 
               outlier.shape=16,
               outlier.size=1,
               notch=TRUE) +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position = 'none') +
  theme(text=element_text(size=10, family="Palatino"))

# Beta convergence scatterplot 
stch36Growth %>% 
  ggplot(aes(x = ln_x, y = growth1)) +
  geom_point(size = 4) + 
  geom_smooth(size = 2) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of TFP in 2003",
       y = "TFP Growth")

4.2 Mobility Analysis

# Select and rename variables
stch36 <- stch36 %>% 
  select(id, ln_r_x, ln_r_y) %>% 
  rename(initial = ln_r_x, final = ln_r_y)

# Identify forward and backward mobility
stch36 <- stch36 %>% 
  mutate(mobility0 = as.factor(final > initial))

# Identify mobility across quadrants
stch36 <- stch36 %>% 
  mutate(mobility = as.factor(case_when(final > 0 & initial > 0 ~ 1, 
                              final > 0 & initial < 0 ~ 2,
                              final < 0 & initial < 0 ~ 3,
                              final < 0 & initial > 0 ~ 4,
                              TRUE ~ 0))
                              )
# Identify observation in the first quadrant (uppper right corner)
stch36 <- stch36 %>% 
  mutate(mobility1 = as.factor(final > 0 & initial > 0))
  
# Identify observation in the second quadrant (lower right corner)
stch36 <- stch36 %>% 
  mutate(mobility2 = as.factor(final > 0 & initial < 0))
# Identify observation in the third quadrant (lower left corner)
stch36 <- stch36 %>% 
  mutate(mobility3 = as.factor(final < 0 & initial < 0))
# Identify observation in the third quadrant (upper left corner)
stch36 <- stch36 %>% 
  mutate(mobility4 = as.factor(final < 0 & initial > 0))

# Print table
stch36

# Create scatterplot 
stch36 %>% 
ggplot(aes(x = final , y = initial, color = mobility0, shape = mobility)) +
  geom_point(size = 4) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Create scatterplot with bivarate density
stch36 %>% 
ggplot(aes(x = final , y = initial)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_gradientn(colors = c("#FFEDA0", "#FEB24C", "#F03B20"))+
  geom_point(size = 2) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Descriptive statistics for continuous variables
kable(descr(stch36))

	initial	final
Mean	-0.21	-0.32
Std.Dev	0.65	0.73
Min	-1.78	-1.70
Median	-0.22	-0.41
Max	1.45	1.96
MAD	0.72	0.67
IQR	0.93	0.90
CV	-0.32	-0.44
Skewness	0.04	0.82
SE.Skewness	0.30	0.30
Kurtosis	-0.14	0.56
N.Valid	62.00	62.00
Pct.Valid	100.00	100.00

# Frequency table for forward and backward mobility
kable(freq(stch36$mobility0))

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
FALSE	38	61	61	61	61
TRUE	24	39	100	39	100
	0	NA	NA	0	100
Total	62	100	100	100	100

# Frequency table for mobility across quadrants
freq(stch36$mobility)

Frequencies   
mobility     
Data frame: stch36   
Type: Factor (unordered)   

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          1      8     12.90          12.90     12.90          12.90
          2      9     14.52          27.42     14.52          27.42
          3     33     53.23          80.65     53.23          80.65
          4     12     19.35         100.00     19.35         100.00
       <NA>      0                               0.00         100.00
      Total     62    100.00         100.00    100.00         100.00

# Simple mobility matrix
kable(with(stch36, ctable(mobility0, mobility, prop = "r")) )

	1	2	3	4	Total
FALSE	6	0	20	12	38
TRUE	2	9	13	0	24
Total	8	9	33	12	62

	1	2	3	4	Total
FALSE	0.16	0.00	0.53	0.32	1
TRUE	0.08	0.38	0.54	0.00	1
Total	0.13	0.15	0.53	0.19	1

5 Sector ISIC36LC-lp

5.1 Sigma and Beta Convergence Analysis

# Data import
lp36 <- read_csv(" ISIC36LC-lp/inputs/dataMobility.csv")
lp36 <- lp36 %>% 
  mutate(id = as.factor(id))

# Compute average annual growth rates (use approximations)
lp36 <- lp36 %>% 
  mutate(growth1 = ((ln_y - ln_x)/6)*100,
         growth2 =  ((y/x)^(1/6)-1)*100)
lp36

# Descriptive statistics of continous variables
lp36Growth <- lp36 %>% 
  select(x, y, ln_x, ln_y, growth1, growth2, r_x, r_y, ln_r_x, ln_r_y) 
kable(descr(lp36Growth))

	x	y	ln_x	ln_y	growth1	growth2	r_x	r_y	ln_r_x	ln_r_y
Mean	42.8	47.9	3.43	3.51	1.48	2.46	1.00	1.00	-0.33	-0.35
Std.Dev	38.1	45.7	0.82	0.85	13.72	14.65	0.89	0.95	0.82	0.85
Min	6.2	5.3	1.82	1.67	-29.71	-25.70	0.14	0.11	-1.94	-2.19
Median	28.5	35.0	3.35	3.56	-0.25	-0.25	0.67	0.73	-0.41	-0.31
Max	196.4	254.5	5.28	5.54	36.64	44.25	4.59	5.31	1.52	1.67
MAD	21.1	26.6	0.88	0.88	10.61	10.57	0.49	0.55	0.88	0.88
IQR	40.8	38.8	1.26	1.14	14.42	14.46	0.95	0.81	1.26	1.14
CV	1.1	1.1	4.20	4.14	0.11	0.17	1.12	1.05	-0.41	-0.42
Skewness	1.7	2.2	0.21	0.09	0.51	0.93	1.71	2.19	0.21	0.09
SE.Skewness	0.3	0.3	0.30	0.30	0.30	0.30	0.30	0.30	0.30	0.30
Kurtosis	3.0	5.7	-0.79	-0.45	0.36	0.79	3.04	5.74	-0.79	-0.45
N.Valid	62.0	62.0	62.00	62.00	62.00	62.00	62.00	62.00	62.00	62.00
Pct.Valid	100.0	100.0	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

# Create long panel and recode year
lp36Panel <- lp36 %>%
  gather(year, ln_tfp, ln_x:ln_y)

# Recode year variable
lp36Panel$year <- fct_recode(lp36Panel$year,
                               "2003" = "ln_x",
                               "2009" = "ln_y")

# Plot basic Boxplot
lp36Panel %>% 
  ggplot(aes(x = year, y = ln_tfp, fill = year)) +
  geom_boxplot(outlier.colour="black", 
               outlier.shape=16,
               outlier.size=1,
               notch=TRUE) +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position = 'none') +
  theme(text=element_text(size=10, family="Palatino"))

# Beta convergence scatterplot 
lp36Growth %>% 
  ggplot(aes(x = ln_x, y = growth1)) +
  geom_point(size = 4) + 
  geom_smooth(size = 2) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of TFP in 2003",
       y = "TFP Growth")

5.2 Mobility Analysis

# Select and rename variables
lp36 <- lp36 %>% 
  select(id, ln_r_x, ln_r_y) %>% 
  rename(initial = ln_r_x, final = ln_r_y)

# Identify forward and backward mobility
lp36 <- lp36 %>% 
  mutate(mobility0 = as.factor(final > initial))

# Identify mobility across quadrants
lp36 <- lp36 %>% 
  mutate(mobility = as.factor(case_when(final > 0 & initial > 0 ~ 1, 
                              final > 0 & initial < 0 ~ 2,
                              final < 0 & initial < 0 ~ 3,
                              final < 0 & initial > 0 ~ 4,
                              TRUE ~ 0))
                              )
# Identify observation in the first quadrant (uppper right corner)
lp36 <- lp36 %>% 
  mutate(mobility1 = as.factor(final > 0 & initial > 0))
  
# Identify observation in the second quadrant (lower right corner)
lp36 <- lp36 %>% 
  mutate(mobility2 = as.factor(final > 0 & initial < 0))
# Identify observation in the third quadrant (lower left corner)
lp36 <- lp36 %>% 
  mutate(mobility3 = as.factor(final < 0 & initial < 0))
# Identify observation in the third quadrant (upper left corner)
lp36 <- lp36 %>% 
  mutate(mobility4 = as.factor(final < 0 & initial > 0))

# Print table
lp36

# Create scatterplot 
lp36 %>% 
ggplot(aes(x = final , y = initial, color = mobility0, shape = mobility)) +
  geom_point(size = 4) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Create scatterplot with bivarate density
lp36 %>% 
ggplot(aes(x = final , y = initial)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_gradientn(colors = c("#FFEDA0", "#FEB24C", "#F03B20"))+
  geom_point(size = 2) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Descriptive statistics for continuous variables
kable(descr(lp36))

	initial	final
Mean	-0.33	-0.35
Std.Dev	0.82	0.85
Min	-1.94	-2.19
Median	-0.41	-0.31
Max	1.52	1.67
MAD	0.88	0.88
IQR	1.26	1.14
CV	-0.41	-0.42
Skewness	0.21	0.09
SE.Skewness	0.30	0.30
Kurtosis	-0.79	-0.45
N.Valid	62.00	62.00
Pct.Valid	100.00	100.00

# Frequency table for forward and backward mobility
kable(freq(lp36$mobility0))

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
FALSE	35	56	56	56	56
TRUE	27	44	100	44	100
	0	NA	NA	0	100
Total	62	100	100	100	100

# Frequency table for mobility across quadrants
freq(lp36$mobility)

Frequencies   
mobility     
Data frame: lp36   
Type: Factor (unordered)   

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          1     11     17.74          17.74     17.74          17.74
          2     11     17.74          35.48     17.74          35.48
          3     31     50.00          85.48     50.00          85.48
          4      9     14.52         100.00     14.52         100.00
       <NA>      0                               0.00         100.00
      Total     62    100.00         100.00    100.00         100.00

# Simple mobility matrix
kable(with(lp36, ctable(mobility0, mobility, prop = "r")) )

	1	2	3	4	Total
FALSE	7	0	19	9	35
TRUE	4	11	12	0	27
Total	11	11	31	9	62

	1	2	3	4	Total
FALSE	0.20	0.00	0.54	0.26	1
TRUE	0.15	0.41	0.44	0.00	1
Total	0.18	0.18	0.50	0.15	1

6 Sector ISIC36LC-acf

6.1 Sigma and Beta Convergence Analysis

# Data import
acf36 <- read_csv(" ISIC36LC-acf/inputs/dataMobility.csv")
acf36 <- acf36 %>% 
  mutate(id = as.factor(id))

# Compute average annual growth rates (use approximations)
acf36 <- acf36 %>% 
  mutate(growth1 = ((ln_y - ln_x)/6)*100,
         growth2 =  ((y/x)^(1/6)-1)*100)
acf36

# Descriptive statistics of continous variables
acf36Growth <- acf36 %>% 
  select(x, y, ln_x, ln_y, growth1, growth2, r_x, r_y, ln_r_x, ln_r_y) 
kable(descr(acf36Growth))

	x	y	ln_x	ln_y	growth1	growth2	r_x	r_y	ln_r_x	ln_r_y
Mean	0.49	1.43	-1.36	-0.81	9.17	12.95	1.00	1.00	-0.65	-1.17
Std.Dev	0.86	6.32	1.02	1.04	23.60	31.48	1.74	4.41	1.02	1.04
Min	0.05	0.10	-3.09	-2.34	-35.61	-29.96	0.09	0.07	-2.38	-2.70
Median	0.23	0.38	-1.45	-0.96	6.70	6.92	0.48	0.27	-0.74	-1.32
Max	5.46	49.88	1.70	3.91	96.25	161.83	11.09	34.82	2.41	3.55
MAD	0.19	0.24	0.92	0.77	21.21	21.85	0.40	0.17	0.92	0.77
IQR	0.30	0.42	1.23	1.02	27.84	29.99	0.62	0.29	1.23	1.02
CV	0.57	0.23	-1.33	-0.78	0.39	0.41	0.57	0.23	-0.64	-1.13
Skewness	3.99	7.28	0.79	1.79	1.10	2.29	3.99	7.28	0.79	1.79
SE.Skewness	0.30	0.30	0.30	0.30	0.30	0.30	0.30	0.30	0.30	0.30
Kurtosis	17.87	52.83	0.43	5.48	2.12	7.42	17.87	52.83	0.43	5.48
N.Valid	62.00	62.00	62.00	62.00	62.00	62.00	62.00	62.00	62.00	62.00
Pct.Valid	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00	100.00

# Create long panel and recode year
acf36Panel <- acf36 %>%
  gather(year, ln_tfp, ln_x:ln_y)

# Recode year variable
acf36Panel$year <- fct_recode(acf36Panel$year,
                               "2003" = "ln_x",
                               "2009" = "ln_y")

# Plot basic Boxplot
acf36Panel %>% 
  ggplot(aes(x = year, y = ln_tfp, fill = year)) +
  geom_boxplot(outlier.colour="black", 
               outlier.shape=16,
               outlier.size=1,
               notch=TRUE) +
  scale_fill_brewer(palette="Blues") +
  theme_classic() +
  theme(legend.position = 'none') +
  theme(text=element_text(size=10, family="Palatino"))

# Beta convergence scatterplot 
acf36Growth %>% 
  ggplot(aes(x = ln_x, y = growth1)) +
  geom_point(size = 4) + 
  geom_smooth(size = 2) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of TFP in 2003",
       y = "TFP Growth")

6.2 Mobility Analysis

# Select and rename variables
acf36 <- acf36 %>% 
  select(id, ln_r_x, ln_r_y) %>% 
  rename(initial = ln_r_x, final = ln_r_y)

# Identify forward and backward mobility
acf36 <- acf36 %>% 
  mutate(mobility0 = as.factor(final > initial))

# Identify mobility across quadrants
acf36 <- acf36 %>% 
  mutate(mobility = as.factor(case_when(final > 0 & initial > 0 ~ 1, 
                              final > 0 & initial < 0 ~ 2,
                              final < 0 & initial < 0 ~ 3,
                              final < 0 & initial > 0 ~ 4,
                              TRUE ~ 0))
                              )
# Identify observation in the first quadrant (uppper right corner)
acf36 <- acf36 %>% 
  mutate(mobility1 = as.factor(final > 0 & initial > 0))
  
# Identify observation in the second quadrant (lower right corner)
acf36 <- acf36 %>% 
  mutate(mobility2 = as.factor(final > 0 & initial < 0))
# Identify observation in the third quadrant (lower left corner)
acf36 <- acf36 %>% 
  mutate(mobility3 = as.factor(final < 0 & initial < 0))
# Identify observation in the third quadrant (upper left corner)
acf36 <- acf36 %>% 
  mutate(mobility4 = as.factor(final < 0 & initial > 0))

# Print table
acf36

# Create scatterplot 
acf36 %>% 
ggplot(aes(x = final , y = initial, color = mobility0, shape = mobility)) +
  geom_point(size = 4) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Create scatterplot with bivarate density
acf36 %>% 
ggplot(aes(x = final , y = initial)) +
  stat_density_2d(aes(fill = ..level..), geom = "polygon") +
  scale_fill_gradientn(colors = c("#FFEDA0", "#FEB24C", "#F03B20"))+
  geom_point(size = 2) + 
  # Add 45 degree line and reference lines
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_hline(yintercept=0, linetype="dashed") + 
  geom_vline(xintercept=0, linetype="dashed") +
  # Add labels that do not overlap
  # geom_text_repel(aes(label = id), size = 3) +
  # Change axis limits
  xlim(-3, 3) + 
  ylim(-3, 3) +
  # Change the theme to viridis colors and black and white background
  scale_color_viridis(discrete = TRUE) +
  theme_bw() +
  # Change the font of the theme
  theme(text=element_text(size=20, family="Palatino")) +
  # Hide the lenged
  theme(legend.position = 'none') +
  # Change plot titles
  labs(x = "Log of Relative TFP in 2009",
       y = "Log of Relative TFP in 2003")

# Descriptive statistics for continuous variables
kable(descr(acf36))

	initial	final
Mean	-0.65	-1.17
Std.Dev	1.02	1.04
Min	-2.38	-2.70
Median	-0.74	-1.32
Max	2.41	3.55
MAD	0.92	0.77
IQR	1.23	1.02
CV	-0.64	-1.13
Skewness	0.79	1.79
SE.Skewness	0.30	0.30
Kurtosis	0.43	5.48
N.Valid	62.00	62.00
Pct.Valid	100.00	100.00

# Frequency table for forward and backward mobility
kable(freq(acf36$mobility0))

	Freq	% Valid	% Valid Cum.	% Total	% Total Cum.
FALSE	44	71	71	71	71
TRUE	18	29	100	29	100
	0	NA	NA	0	100
Total	62	100	100	100	100

# Frequency table for mobility across quadrants
freq(acf36$mobility)

Frequencies   
mobility     
Data frame: acf36   
Type: Factor (unordered)   

              Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
----------- ------ --------- -------------- --------- --------------
          1      1      1.61           1.61      1.61           1.61
          2      4      6.45           8.06      6.45           8.06
          3     45     72.58          80.65     72.58          80.65
          4     12     19.35         100.00     19.35         100.00
       <NA>      0                               0.00         100.00
      Total     62    100.00         100.00    100.00         100.00

# Simple mobility matrix
kable(with(acf36, ctable(mobility0, mobility, prop = "r")) )

	1	2	3	4	Total
FALSE	1	0	31	12	44
TRUE	0	4	14	0	18
Total	1	4	45	12	62

	1	2	3	4	Total
FALSE	0.02	0.00	0.70	0.27	1
TRUE	0.00	0.22	0.78	0.00	1
Total	0.02	0.06	0.73	0.19	1