Analiza determinantów zanieczyszczeń powietrza
Jakub Łukasiewicz
6 03 2021
wstęp
Zagadnienia ochrony środowiska stają się coraz bardziej popularne w kulturze oraz mediach. Wraz ze zwiększoną ilościa informacji na ten temat rośnie również świadomość społeczeństwa. Znieczyszczenie środowiska gwałtownie wzrasta wraz ze wzrostem konsumpcjonizmu oraz światowej liczby ludności. W 2016 roku Unia Europejska wuprodukowała blisko 2,5 miliarda ton odpadów , a średni europejczyk produkuje blisko 500 kilogramów odpoadów rocznie (stan na 2018r.). Równie istotnym zagadnieniem jest również zanieczyszczenie powietrza oraz emisja CO2. Wysoko rozwinięte państwa, świadome istniejącego zagrożenia starają się ograniczyć emisje zanieczyszeń do środowiska. Coraz więcej państwo decyduje się na wprowadzenie długoterminowych planów mających na celu osiągnięcie neutralności klimatycznej. W tym celu maksymalnie ogranicza się emisję gazów cieplarnianych we wsystkich możliwych obszarach np. w przemyśle, transporcie itd.. Dodatkowo stawiając na sadzenie drzew, które pochłaniają część gazów np. CO2. Unia Europejska założyła, że osiągnie ten plan do 2050 roku. Ważnym elementem tego planu jest m.in. ograniczanie udziału elektrowni konwencjonalnych w produkcji energii elektrycznej. Coraz więcej stawia się na źródła odnawialne np. farmy wiatrowe oraz panele fotowoltaiczne. Dla przykładu, elektrowania PAK KWB Konin zamierza do 2030 zamkąć elektrownie węglowe i zastąpić je w pełni odnawialnymi źródłami. W mojej pracy chciałbym skupić się na zbadanie determinantów emisji zanieczyszczeń powietrza na przykładzie 14 regionów administracyjnych Czech. W mojej analizie użyję przestrzennego modelu ekonometrycznego badając dyfujzę pomiędzy regionami. Przenalizuję wpływ transportu, zużycia energi oraz PKB per capita na emisję. Dodatkowyn elementem, który wprowadzę są wydatki na badania i rozwój (R&D). Poniżej mapa Czech z podziałem na regiony.
library(spdep)## Loading required package: sp## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Loading required package: sf## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1library(rgdal)## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: C:/Program Files/R/R-4.0.4/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Program Files/R/R-4.0.4/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## Overwritten PROJ_LIB was C:/Program Files/R/R-4.0.4/library/rgdal/projlibrary(maptools)## Checking rgeos availability: TRUElibrary(sp)
library(RColorBrewer)
library(classInt)
library(GISTools)## Loading required package: MASS## Loading required package: rgeos## rgeos version: 0.5-5, (SVN revision 640)
## GEOS runtime version: 3.8.0-CAPI-1.13.1
## Linking to sp version: 1.4-5
## Polygon checking: TRUElibrary(maps)##
## Attaching package: 'maps'## The following object is masked from 'package:GISTools':
##
## map.scale#Wczytanie mapy Czech
CZE <- readOGR(".","CZE_adm1")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Kuba\Documents\ER\CZE", layer: "CZE_adm1"
## with 14 features
## It has 9 fields
## Integer64 fields read as strings: ID_0 ID_1plot(CZE)
#Wczytanie danych
enr <- read.csv("energy.csv", sep = ";")
stat <- read.csv("data.csv", sep = ";")
dane <- merge(x=stat, y=enr, by.x = 'ID_1', by.y = 'ID')
View(dane)
dane$R.D <- dane$R.D*0.045*1000000
dane$R.Dpc <- dane$R.D/dane$PopulationPRZEGLĄD LITERATURY
Inspiracją do napisania tego artykułu była praca Poon, Casas and He “The Impact of Energy, Transport, and Trade on Air Pollution in China”. Autorzy tej pracy zbadali wpływ zużycia energii, transportu oraz handlu na zanieczyszczenie powietrza. Badanie zostało przprowadzone dla prowincji chińskich dla 2004 roku. W artykule wykzano teorię Kuznetsova o odwrotnej krzywej U, potwierdza się dla SO2. Ilość konsumowanej energii elektrycznej oraz ilośc samochodów zwiększa poziom znaieczyszczeń. Odwrotny wpływ wiąże ze otwartość handlowa. Liu i inni (2017) przprowadzili podobną analizę również dla regionów administracyjnych Chin. Zbadnao jednak znacznie większą ilość zmiennych. Dodano zmienne demograficzne takie jka populacja czy też urbanizacja. Podobnie jak w poprzedniej pracy do analizy wykorzystano również transport w postaci gęstości samochodów. W pracy potwierdzono występowanie efektu dyfuzji badanego zjawiska. Wskazano również na dodatni wpływ urbanizacji na zanieczyszczenie powietrza. W opozycji do innych tekstów w tej tematyce, nie wykazano istotnego wpływu PKB per capita na ilość zanieczyszczeń. W pracy Wu i inni (2019) do badania nad zanieczyszczeniem powietrza użyto modeli Durbina. W tym artykule zbadany został wpływ regulacji na emisję zanieczyszczeń. Zgodnie z wynikami artykuły, wykazany został pozytywny wpływ regulacji na emisję zanieczyszczeń. Do kontroli zmiennych w tej pracy zostały wykorzystane m.in. konsumpcja energii elektrycznej oraz PKB per capita.
WSTĘP TEORETYCZNY
Jednym z czynników kształtujących emisję znaieczyszczeń powietrza jest wzrost gospodarczy. Zgodnie z teorią Kuznetsa wraz ze wzrostem gospodarczym, wzrasta również poziom zanieczyszczeń. Krzywa Kuznetsa ma kształt odwróconej litery U, co implikuje, że po osiągnięciu pewnego poziomu gospodarczego poziom emitowanych zanieczyszczeń spada (Poon, 2006). Jest to zbieżne z inwestcjami w badania i rozwój, ponieważ rozwój technologiczny wiąże się często z bardziej oszczędnymi i ekologicznymi metodami produkcji.
Konsumpcja energii jest ściśle związana z produkcją i przemysłem. Duże zużycie może oznaczć istnienie dużych zakładów wytwórczych itp. w regionie. Duża konsumpcja wiąże się również z produkcją energii, która często jest źródłem dużej ilości zanieczyszczeń.
Mniej powszechny wpływ na zanieczyszczenie powietrza przywiązuje się do transportu. Problem ten jest jednak dość dobrze obserwowany w dużych miastach. Jak wskazuje Poon w sowjej pracy, może być to jeden z czynników wpływających na znieczyszczenie. BAZA DANYCH ORAZ METODOLOGIA
Zbiór danych zawiera 14 regionów administracyjnych Republiki Czeskiej za rok 2018. Z uwagi na brak dostępnych danych nie było możliwości przejścia na bardziej szczegółowy model pod względem podziału administracyjnego. Wszytskie zmienne zostały przedstawione w tabeli poniżej.meta <- read.csv("meta.csv", sep = ";")
meta## Zmienna Opis
## 1 CO2 Emisja dwutlenku węgla dla regionu w tonach rocznie.
## 2 SO2 Emisja dwutleneku siarki dla regionu w tonach rocznie.
## 3 PKBpc Nominalny Prodkut Krajowy Brutto per capita w dolarach
## 4 CAR Liczba zarejestrownaych samochodów osobowych
## 5 R.D Wydatki na badania i rozwój w USD
## 6 ENERGY Roczne zużycie energii w MWh
## Źródło
## 1 Czeski Urząd Statystyczny
## 2 Czeski Urząd Statystyczny
## 3 Czeski Urząd Statystyczny
## 4 Ministerstwo Transportu Czech
## 5 Czeski Urząd Statystyczny
## 6 Czeski Urząd Regulacji EnergiiNiektóre zmienne w analizie zostały zlogarytmowane, co oznaczone zostało przez dopisek log. W poniższej analiznie zdecydowałem się na model przestrzenny ze względu na Możliwą dyfuzję zanieczyszczeń i ich przepływ z regionów mocno zanieczyszczonych do ich bardziej ekologicznych sąsiadów.
ANALIZA ORAZ REZULTATY
Analizę przestrzenną danego zjawiska rozpocząłem od przedstawienia na mapie danych dotyczących emisji SO2 oraz konsumpcji energii w 2018 roku. Jak mozęmy zauważyć, regiony znacznie różnią się pod tym względem między sobą. Widoczne są pewne zależności pomiędzy emisją SO2 a konsumpcją energii.
#wykres SO2
zmienna <- dane$SO2
zmienna## [1] 25638.5 3242.0 1278.1 6527.8 3417.3 1729.6 1027.6 15755.6 3177.4
## [10] 5897.3 3876.7 258.5 15181.5 3173.8summary(zmienna)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 258.5 2090.7 3329.7 6441.6 6370.2 25638.5mean(zmienna)## [1] 6441.55przedziały<-6
kolory<-brewer.pal(przedziały, "BuPu") # wybór kolorów
klasy<-classIntervals(zmienna, przedziały, style="fixed",
fixedBreaks=c(5000, 10000, 15000,20000,25000,30000))## Warning in classIntervals(zmienna, przedziały, style = "fixed", fixedBreaks =
## c(5000, : variable range greater than fixedBreakstabela.kolorów<-findColours(klasy, kolory)
plot(CZE, col=tabela.kolorów)
plot(CZE, lwd=2, add=TRUE)
legend("bottomleft", legend=names(attr(tabela.kolorów,
"table")), fill=attr(tabela.kolorów, "palette"), cex=1, bty="n")
title(main="Emisja SO2 w tonach w 2018 r.")
#Wykres dla konsumpcji ee
zmienna <- dane$TOTAL
zmienna## [1] 6134820 3154770 5434045 1535949 3601581 2747153 2549863 8189796 3213656
## [10] 2488105 3080488 6057999 8036205 3287290summary(zmienna)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1535949 2830486 3250473 4250837 5902010 8189796mean(zmienna)## [1] 4250837przedziały<-8
kolory<-brewer.pal(przedziały, "BuPu") # wybór kolorów
klasy<-classIntervals(zmienna, przedziały, style="fixed",
fixedBreaks=c(2000000, 3000000, 4000000, 5000000, 6000000, 7000000, 8000000, 9000000))## Warning in classIntervals(zmienna, przedziały, style = "fixed", fixedBreaks =
## c(2e+06, : variable range greater than fixedBreakstabela.kolorów<-findColours(klasy, kolory)
plot(CZE, col=tabela.kolorów)
plot(CZE, lwd=2, add=TRUE)
legend("bottomleft", legend=names(attr(tabela.kolorów,
"table")), fill=attr(tabela.kolorów, "palette"), cex=1, bty="n")
title(main="Konsumpcja energii elektrycznej w MWh w 2018 r.")
W celu zbadania autokorelacji przestrzennej dla zmiennej SO2 policzono statystykę I Morana, która sugeruje nam ujemną relację,
cont.nb<-poly2nb(as(CZE, "SpatialPolygons"))
cont.listw<-nb2listw(cont.nb, style="W")
zmienna <- dane$SO2
wynik02<-moran.test(zmienna, cont.listw)
wynik03<-moran.test(zmienna, cont.listw, randomisation=FALSE)
zmienna<-as.data.frame(scale(dane$SO2))# standaryzacja
moran.plot(zmienna$V1, cont.listw, pch=19, labels=as.character(dane$Region)) 
Na poniższym wykresie przedstawiono przynależność regionów do ćwiartek punktowej statysyki Morana.
#przygotowanie danych
x<-dane$SO2
zx<-scale(x)
wzx<-lag.listw(cont.listw, zx)
cond1<-ifelse(zx>=0 & wzx>=0, 1,0) # I cwiartka
cond2<-ifelse(zx>=0 & wzx<0, 2,0) # II cwiartka
cond3<-ifelse(zx<0 & wzx<0, 3,0) # III cwiartka
cond4<-ifelse(zx<0 & wzx>=0, 4,0) # IV cwiartka
cond.all<-cond1+cond2+cond3+cond4
cond.all## [,1]
## [1,] 1
## [2,] 3
## [3,] 3
## [4,] 1
## [5,] 4
## [6,] 3
## [7,] 4
## [8,] 2
## [9,] 4
## [10,] 3
## [11,] 4
## [12,] 4
## [13,] 2
## [14,] 4cond<-as.data.frame(cond.all)
is.data.frame(cond)## [1] TRUE# wykres
brks<-c(1,2,3,4)
cols<-c("grey25", "grey60", "grey85", "grey45")
plot(CZE, col=cols[findInterval(cond$V1, brks)])
legend("bottomleft", legend=c("I ćw - HH - wysokie otoczone wysokimi", "II ćw - LH - niskie otoczone wysokimi", "III ćw - LL - niskie otoczone niskimi", "IV ćw - HL - wysokie otoczone niskimi"), fill=cols, bty="n", cex=0.85)
title(main="Przynależność regionów do ćwiartek
z wykresu punktowego Morana")
# Model
Przeprowadzony test Morana wskazuje nam, że jest potrzeba zastosowania modelu przestrzennego.
#Przygotowanie danych
sub<-dane
sub$y<-sub$SO2
sub$x1<-sub$CARS
sub$x2<-sub$GDPpc
sub$x3<-sub$R.Dpc
sub$x4<-sub$TOTAL
sub$x5<-sub$R.D
# model liniowy
model.lm<-lm(y~x1+x2+x3+x4, data=sub)
summary(model.lm)##
## Call:
## lm(formula = y ~ x1 + x2 + x3 + x4, data = sub)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4865.2 -2755.9 331.2 1857.0 7999.5
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.050e+03 5.063e+03 -0.800 0.4444
## x1 1.879e-02 3.102e-02 0.606 0.5597
## x2 2.658e-01 3.506e-01 0.758 0.4677
## x3 -3.075e+01 1.497e+01 -2.055 0.0701 .
## x4 1.855e-03 1.974e-03 0.940 0.3719
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4240 on 9 degrees of freedom
## Multiple R-squared: 0.7674, Adjusted R-squared: 0.664
## F-statistic: 7.424 on 4 and 9 DF, p-value: 0.006286lm.morantest(model.lm, cont.listw)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = y ~ x1 + x2 + x3 + x4, data = sub)
## weights: cont.listw
##
## Moran I statistic standard deviate = 1.2561, p-value = 0.1045
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.13396504 -0.08217417 0.02960671Do zbadania zależności przestrzennej zostanie wykorzystany prosty model jednokomponentowy opóźnienia przestrzennego. Zmienną objaśnianą jest zmienna SO2, czyli emisja dwutlenku siarki w tonach. Do objaśnienia danego zjawiska wybrano ilość zarejestrowanych aut, PKB per capita, R&D per capita oraz konsumpcję energii. Zbudowane zostały dwa modele, jeden z noaminalnymi wartościami oraz drugi z wartościami częsciowo zlogarytmowanymi. W pierwszym modelu jedynie wydatki na badania per capita okazały się istotną zmienna na poziomie istotności 10%, nieznacznie powyżej 10% poziomu istotności znalazł się współczynnik rho. Zmienna R&D w tym modelu ma ujemny wpływ na wielkość emisji co pokrywa się z przewidywaniami. W drugim modelu niestotne okazała się zmienna dotycząca liczby zarejestrowanych aut oraz współczynnik opóźnienia rho. Kierunki zależnoci zmiennych zdają się potwierdzać przedstawioną wcześniej teorię. Zarówno ilość aut, jak i wielkość zużycia energii zwiększa emisję dwutlenku siarki do atmosfery. Tak jak w poprzednim modelu, również i tutaj badania i rozwój mają pozytywny wpływ na redukcję zanieczyszczeń. Ujemny współczynnik przy PKB per capita oznacza, że Czechy należą do grona państw rozwiniętych, w których to wzrost gospodarczy nie przyczynia się do większej emisji gazów do atmosfery.
model.lag<-lagsarlm(y~x1+x2+x3+x4, data=sub, cont.listw, tol.solve=3e-30)
wynik.lag<-summary(model.lag)
wynik.lag##
## Call:lagsarlm(formula = y ~ x1 + x2 + x3 + x4, data = sub, listw = cont.listw,
## tol.solve = 3e-30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3960.428 -1831.865 48.837 710.162 7371.501
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.3936e+03 3.9253e+03 -1.6288 0.10335
## x1 2.3586e-02 2.2706e-02 1.0387 0.29892
## x2 1.3656e-01 2.6231e-01 0.5206 0.60266
## x3 -3.1850e+01 1.0924e+01 -2.9156 0.00355
## x4 1.9603e-03 1.4318e-03 1.3691 0.17098
##
## Rho: 0.36912, LR test value: 2.2708, p-value: 0.13183
## Asymptotic standard error: 0.22898
## z-value: 1.6121, p-value: 0.10695
## Wald statistic: 2.5987, p-value: 0.10695
##
## Log likelihood: -132.5697 for lag model
## ML residual variance (sigma squared): 9450200, (sigma: 3074.1)
## Number of observations: 14
## Number of parameters estimated: 7
## AIC: 279.14, (AIC for lm: 279.41)
## LM test for residual autocorrelation
## test value: 0.9907, p-value: 0.31957bptest.sarlm(model.lag) # heteroskedasticity##
## studentized Breusch-Pagan test
##
## data:
## BP = 5.7451, df = 4, p-value = 0.219model.lag<-lagsarlm(y~log(x1)+x2+log(x3)+log(x4), data=sub, cont.listw, tol.solve=3e-30)
wynik.lag<-summary(model.lag)
wynik.lag##
## Call:lagsarlm(formula = y ~ log(x1) + x2 + log(x3) + log(x4), data = sub,
## listw = cont.listw, tol.solve = 3e-30)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4163.47 -1434.99 263.87 570.95 6495.97
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.8880e+05 3.0613e+04 -6.1673 6.946e-10
## log(x1) 4.7376e+03 6.8119e+03 0.6955 0.48675
## x2 -2.7532e-01 1.6208e-01 -1.6987 0.08938
## log(x3) -6.2585e+03 1.3136e+03 -4.7645 1.894e-06
## log(x4) 1.1378e+04 5.2026e+03 2.1869 0.02875
##
## Rho: 0.2985, LR test value: 1.5885, p-value: 0.20754
## Asymptotic standard error: 0.21241
## z-value: 1.4053, p-value: 0.15993
## Wald statistic: 1.9749, p-value: 0.15993
##
## Log likelihood: -130.3673 for lag model
## ML residual variance (sigma squared): 6998400, (sigma: 2645.4)
## Number of observations: 14
## Number of parameters estimated: 7
## AIC: 274.73, (AIC for lm: 274.32)
## LM test for residual autocorrelation
## test value: 1.8356, p-value: 0.17547bptest.sarlm(model.lag) # heteroskedasticity##
## studentized Breusch-Pagan test
##
## data:
## BP = 1.3413, df = 4, p-value = 0.8543PODSUMOWANIE
Przedsatwiona analiza potwierdza przedstawione w dziale teoretycznym założenia oraz intiucję. Zarówno rosnąca ilość aut, jak i zwiększona konumpcja energii negatywnie wpływa na środowisko i przyczynia się do zwiększenia emisji gazów cieplarnianych do środowiska. Pozytywny wpływ na emisję, prowadzący do jej ograniczenia ma zaproponowana przeze mnie zmienna dotycząca wydatków na badnia i rozwój (R&D). Oznacza to, że środki te mogą przczyniać się do powstawaniach nowych technologi, które okazują się nie tylko bardziej produktywne ale również mniej szkodliwe dla społeczeństwa. Z uwagi na małą liczbę obserwacji, wynik modelu może być niedokładny. Głównym problemem okazuję się mała dostępność danych dot. zużycia energii czy też emisji zanieczyszczeń na niższym poziomie administracyjnym. Analizę należałoby powtórzyć po pojawieniu się odpowiednich danych umożliwiających dokładniejsze zbadanie zjawiska.
BIBLIOGRAFIA
Jessie P. H. Poon, Irene Casas & Canfei He (2006) The Impact of Energy, Transport, and Trade on Air Pollution in China, Eurasian Geography and Economics, 47:5, 568-584, DOI: 10.2747/1538-7216.47.5.568
Liu, H., Fang, C., Zhang, X., Wang, Z., Bao, C., & Li, F. (2017). The effect of natural and anthropogenic factors on haze pollution in Chinese cities: A spatial econometrics approach. Journal of Cleaner Production, 165, 323–333. doi:10.1016/j.jclepro.2017.07.127
Wu, X., Gao, M., Guo, S., & Li, W. (2019). Effects of environmental regulation on air pollution control in China: a spatial Durbin econometric analysis. Journal of Regulatory Economics. doi:10.1007/s11149-019-09384-x
https://300gospodarka.pl/explainer/co-to-jest-neutralnosc-klimatyczna