It’s just playing around with quantile regression on a pooled dataset.
rm(list = ls())
setwd("/Users/vancam/Documents/WAIKATO-Thesis/VESuse/CVan/TFP_OwnGO")
load("/Users/vancam/Documents/WAIKATO-Thesis/VESuse/CVan/Final_forpaper/export_ublfinal.rda")
library(tidyverse)
## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr
## Conflicts with tidy packages ----------------------------------------------
## filter(): dplyr, stats
## lag(): dplyr, stats
export <- filter(export_ublfinal, exportvol >0)
export <- mutate(export, loglabor = log(labor))
export <- mutate(export, size = labor)
export <- mutate(export, sizesq = labor*labor)
export <- mutate(export, logsize = log(size))
attach(export)
## The following object is masked _by_ .GlobalEnv:
##
## export
Y <- cbind(exportvol)
X <- cbind(size, scale, factor(ownership), Hlinkage, industrialzone)
library(quantreg)
## Loading required package: SparseM
##
## Attaching package: 'SparseM'
## The following object is masked from 'package:base':
##
## backsolve
ols <- lm(Y~X, data = export)
summary(ols)
##
## Call:
## lm(formula = Y ~ X, data = export)
##
## Residuals:
## Min 1Q Median 3Q Max
## -102717 -2104 -494 1098 394890
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6567.1362 363.0450 18.089 < 2e-16 ***
## Xsize 11.0497 0.1427 77.412 < 2e-16 ***
## Xscale 243.4551 4.2749 56.950 < 2e-16 ***
## X -4112.6279 309.4852 -13.289 < 2e-16 ***
## XHlinkage -6248.6533 317.6572 -19.671 < 2e-16 ***
## Xindustrialzone 1253.8619 178.6808 7.017 2.34e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10040 on 19161 degrees of freedom
## Multiple R-squared: 0.4292, Adjusted R-squared: 0.429
## F-statistic: 2881 on 5 and 19161 DF, p-value: < 2.2e-16
#Quantiles
quan10 <- rq(Y~X, data = export, tau = 0.1)
quan20 <- rq(Y~X, data = export, tau = 0.2)
quan30 <- rq(Y~X, data = export, tau = 0.3)
quan40 <- rq(Y~X, data = export, tau = 0.4)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
quan50 <- rq(Y~X, data = export, tau = 0.5)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
quan60 <- rq(Y~X, data = export, tau = 0.6)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
quan70 <- rq(Y~X, data = export, tau = 0.7)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
quan80 <- rq(Y~X, data = export, tau = 0.8)
quan90 <- rq(Y~X, data = export, tau = 0.9)
summary(quan10)
## Warning in summary.rq(quan10): 459 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.1, data = export)
##
## tau: [1] 0.1
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 15.24960 4.72373 3.22830 0.00125
## Xsize 0.16916 0.02985 5.66738 0.00000
## Xscale 1.06311 0.44298 2.39988 0.01641
## X -3.58615 4.92986 -0.72743 0.46697
## XHlinkage -10.21246 1.91136 -5.34302 0.00000
## Xindustrialzone 7.12891 3.45440 2.06372 0.03906
summary(quan20)
## Warning in summary.rq(quan20): 1146 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.2, data = export)
##
## tau: [1] 0.2
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 78.06838 6.18869 12.61468 0.00000
## Xsize 0.97188 0.09234 10.52525 0.00000
## Xscale 3.10448 5.54116 0.56026 0.57531
## X -60.00353 7.59780 -7.89749 0.00000
## XHlinkage -38.76316 4.18777 -9.25627 0.00000
## Xindustrialzone 13.09452 5.48949 2.38538 0.01707
summary(quan30)
## Warning in summary.rq(quan30): 1402 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.3, data = export)
##
## tau: [1] 0.3
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 211.55002 23.24330 9.10155 0.00000
## Xsize 2.11458 0.17570 12.03486 0.00000
## Xscale 25.18240 13.86862 1.81578 0.06942
## X -173.05782 25.48325 -6.79104 0.00000
## XHlinkage -110.87141 9.08251 -12.20713 0.00000
## Xindustrialzone 9.45440 2.65725 3.55796 0.00037
summary(quan40)
## Warning in summary.rq(quan40): 986 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.4, data = export)
##
## tau: [1] 0.4
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 400.76055 33.85094 11.83898 0.00000
## Xsize 3.32785 0.27105 12.27753 0.00000
## Xscale 82.63497 31.40032 2.63166 0.00850
## X -333.55664 33.05641 -10.09053 0.00000
## XHlinkage -205.72697 15.88470 -12.95127 0.00000
## Xindustrialzone 10.29999 13.81507 0.74556 0.45594
summary(quan50)
## Warning in summary.rq(quan50): 815 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.5, data = export)
##
## tau: [1] 0.5
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 474.90011 51.86717 9.15608 0.00000
## Xsize 4.34419 0.33993 12.77957 0.00000
## Xscale 231.47953 35.34815 6.54856 0.00000
## X -371.36910 52.19331 -7.11526 0.00000
## XHlinkage -309.53603 13.54900 -22.84567 0.00000
## Xindustrialzone -6.43663 5.61223 -1.14689 0.25144
summary(quan60)
## Warning in summary.rq(quan60): 729 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.6, data = export)
##
## tau: [1] 0.6
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 553.68819 82.60589 6.70277 0.00000
## Xsize 6.79834 0.42361 16.04876 0.00000
## Xscale 284.80673 32.20164 8.84448 0.00000
## X -384.78890 83.04574 -4.63346 0.00000
## XHlinkage -488.90205 24.24131 -20.16813 0.00000
## Xindustrialzone 13.52652 26.99804 0.50102 0.61636
summary(quan70)
## Warning in summary.rq(quan70): 377 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.7, data = export)
##
## tau: [1] 0.7
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 771.85644 85.44348 9.03353 0.00000
## Xsize 9.64132 0.53984 17.85946 0.00000
## Xscale 384.62417 51.91733 7.40840 0.00000
## X -472.17945 80.27040 -5.88236 0.00000
## XHlinkage -754.30673 48.25191 -15.63268 0.00000
## Xindustrialzone 59.36857 45.24502 1.31216 0.18948
summary(quan80)
## Warning in summary.rq(quan80): 290 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.8, data = export)
##
## tau: [1] 0.8
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 1127.12436 164.83562 6.83787 0.00000
## Xsize 14.11493 0.62444 22.60432 0.00000
## Xscale 529.31077 70.39604 7.51904 0.00000
## X -530.19854 160.72485 -3.29880 0.00097
## XHlinkage -1255.08062 93.06621 -13.48589 0.00000
## Xindustrialzone 168.17722 94.16582 1.78597 0.07412
summary(quan90)
## Warning in summary.rq(quan90): 207 non-positive fis
##
## Call: rq(formula = Y ~ X, tau = 0.9, data = export)
##
## tau: [1] 0.9
##
## Coefficients:
## Value Std. Error t value Pr(>|t|)
## (Intercept) 2396.02053 249.58022 9.60020 0.00000
## Xsize 19.86599 1.32495 14.99376 0.00000
## Xscale 853.12093 208.03324 4.10089 0.00004
## X -786.65801 224.64980 -3.50171 0.00046
## XHlinkage -2567.17855 166.23968 -15.44263 0.00000
## Xindustrialzone 391.77853 67.03016 5.84481 0.00000
library(stargazer)
##
## Please cite as:
## Hlavac, Marek (2015). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2. http://CRAN.R-project.org/package=stargazer
stargazer(ols,quan10, quan20, quan30, quan40, quan50, quan60, quan70, quan80, quan90, title = "Export and firms' size - Quantile regression", type = "text")
##
## Export and firms' size - Quantile regression
## ==============================================================================================================================================================
## Dependent variable:
## ------------------------------------------------------------------------------------------------------------------------------------------
## Y
## OLS quantile
## regression
## (1) (2) (3) (4) (5) (6) (7) (8) (9) (10)
## --------------------------------------------------------------------------------------------------------------------------------------------------------------
## Xsize 11.050*** 0.169*** 0.972*** 2.115*** 3.328*** 4.344*** 6.798*** 9.641*** 14.115*** 19.866***
## (0.143) (0.030) (0.092) (0.176) (0.271) (0.340) (0.424) (0.540) (0.624) (1.325)
##
## Xscale 243.455*** 1.063** 3.104 25.182* 82.635*** 231.480*** 284.807*** 384.624*** 529.311*** 853.121***
## (4.275) (0.443) (5.541) (13.869) (31.400) (35.348) (32.202) (51.917) (70.396) (208.033)
##
## X -4,112.628*** -3.586 -60.004*** -173.058*** -333.557*** -371.369*** -384.789*** -472.179*** -530.199*** -786.658***
## (309.485) (4.930) (7.598) (25.483) (33.056) (52.193) (83.046) (80.270) (160.725) (224.650)
##
## XHlinkage -6,248.653*** -10.212*** -38.763*** -110.871*** -205.727*** -309.536*** -488.902*** -754.307*** -1,255.081*** -2,567.179***
## (317.657) (1.911) (4.188) (9.083) (15.885) (13.549) (24.241) (48.252) (93.066) (166.240)
##
## Xindustrialzone 1,253.862*** 7.129** 13.095** 9.454*** 10.300 -6.437 13.527 59.369 168.177* 391.779***
## (178.681) (3.454) (5.489) (2.657) (13.815) (5.612) (26.998) (45.245) (94.166) (67.030)
##
## Constant 6,567.136*** 15.250*** 78.068*** 211.550*** 400.761*** 474.900*** 553.688*** 771.856*** 1,127.124*** 2,396.021***
## (363.045) (4.724) (6.189) (23.243) (33.851) (51.867) (82.606) (85.443) (164.836) (249.580)
##
## --------------------------------------------------------------------------------------------------------------------------------------------------------------
## Observations 19,167 19,167 19,167 19,167 19,167 19,167 19,167 19,167 19,167 19,167
## R2 0.429
## Adjusted R2 0.429
## Residual Std. Error 10,035.630 (df = 19161)
## F Statistic 2,881.167*** (df = 5; 19161)
## ==============================================================================================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## Test for the difference in the intercepts between quantiles
anova(quan10,quan20, joint = FALSE)
## Warning in summary.rq(x, se = se, covariance = TRUE): 459 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 1146 non-positive fis
## Quantile Regression Analysis of Deviance Table
##
## Model: Y ~ X
## Tests of Equality of Distinct Slopes: tau in { 0.1 0.2 }
##
## Df Resid Df F value Pr(>F)
## Xsize 1 38333 105.8018 < 2.2e-16 ***
## Xscale 1 38333 0.1507 0.6979
## X 1 38333 68.8906 < 2.2e-16 ***
## XHlinkage 1 38333 55.8305 8.06e-14 ***
## Xindustrialzone 1 38333 1.7608 0.1845
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(quan10,quan30, joint = FALSE)
## Warning in summary.rq(x, se = se, covariance = TRUE): 459 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 1402 non-positive fis
## Quantile Regression Analysis of Deviance Table
##
## Model: Y ~ X
## Tests of Equality of Distinct Slopes: tau in { 0.1 0.3 }
##
## Df Resid Df F value Pr(>F)
## Xsize 1 38333 138.5613 < 2.2e-16 ***
## Xscale 1 38333 3.1215 0.07727 .
## X 1 38333 51.4197 7.594e-13 ***
## XHlinkage 1 38333 141.1708 < 2.2e-16 ***
## Xindustrialzone 1 38333 0.4265 0.51373
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(quan10,quan20,quan30,quan40,quan50,quan60,quan70,quan80,quan90, joint = FALSE)
## Warning in summary.rq(x, se = se, covariance = TRUE): 459 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 1146 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 1402 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 986 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 815 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 729 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 377 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 290 non-positive fis
## Warning in summary.rq(x, se = se, covariance = TRUE): 207 non-positive fis
## Quantile Regression Analysis of Deviance Table
##
## Model: Y ~ X
## Tests of Equality of Distinct Slopes: tau in { 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 }
##
## Df Resid Df F value Pr(>F)
## Xsize 8 172495 68.4291 < 2.2e-16 ***
## Xscale 8 172495 14.2249 < 2.2e-16 ***
## X 8 172495 22.8420 < 2.2e-16 ***
## XHlinkage 8 172495 88.0873 < 2.2e-16 ***
## Xindustrialzone 8 172495 7.8102 1.523e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Ploting the data
kevin <- rq(Y~X, tau = seq(0.05, 0.95, by = 0.05), data = export)
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique
kevin.plot <- summary(kevin)
## Warning in summary.rq(xi, U = U, ...): 459 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 1190 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 1146 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 1240 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 1402 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 989 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 986 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 749 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 815 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 987 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 729 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 573 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 377 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 282 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 290 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 164 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 207 non-positive fis
## Warning in summary.rq(xi, U = U, ...): 2 non-positive fis
plot(kevin.plot)
phanvi<- rq(exportvol ~ size, data = export, tau=c(0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9))
van <- data.frame(t(coef(phanvi)))
colnames(van) <- c("intercept", "slope")
van$quantiles <- rownames(van)
van
## intercept slope quantiles
## tau= 0.1 7.077566 0.1982763 tau= 0.1
## tau= 0.2 2.019794 0.9982005 tau= 0.2
## tau= 0.3 -6.844961 2.2635659 tau= 0.3
## tau= 0.4 -14.421348 3.9044944 tau= 0.4
## tau= 0.5 -21.714286 6.1428571 tau= 0.5
## tau= 0.6 -34.182082 9.3986217 tau= 0.6
## tau= 0.7 -38.663805 13.8784514 tau= 0.7
## tau= 0.8 34.365471 19.8049327 tau= 0.8
## tau= 0.9 373.484211 33.2315789 tau= 0.9
p<- qplot(x=size, y=exportvol, data=export)
p+ geom_abline(aes(intercept=intercept, slope=slope,colour=quantiles), data=van)