podatki <- read.table("D:/Faks/Ljubljana_Ekonomska/Magisterij/1.letnik/Ekonomske politike EU/Seminarska/end projekt/PANEL/data_use.csv", header = TRUE, sep=";", dec=",")

PODATKI: - Podatki vključujejo 27 evropskih držav. Za vsako državo so pridobljeni podatki za 10 zaporednih let (2013 - 2022)
- Iz podatkov ne bom odstranjeval enot, saj želim ohraniti celovit vzorec vseh 27 enot.

SPREMENLJIVKE:

- REN (odvisna spremenljivka):

Razlaga: Delež energije iz obnovljivih virov predstavlja odstotek končne porabe energije, ki prihaja iz obnovljivih virov. To vključuje električno energijo, toploto in energijo za transport, ki izvirajo iz obnovljivih virov. (% celotne končne porabe energije države).
Napoved: Delež energije iz obnovljivih virov sem izbral kot odvisno spremenljivko, ker me zanima vpliv pojasnjevalnih spremenljivk na stopnjo deleža obnovljive energije.
Vir: Eurostat - https://ec.europa.eu/eurostat/databrowser/view/nrg_ind_ren/default/table?lang=en

- BDP:
Razlaga: BDP na prebivalca (konstantni dolarji ZDA leta 2015).
Napoved: Bruto domači proizvod na prebivalca sem izbral kot glavni ekonomski kazalec gospodarskega razvoja države. Predvidevam, da bo višji BDP na prebivalca vodil v višji delež obnovljive energije, saj ima država več finančnih sredstev za spodbujanje razvoja različnih energetskih virov.
Vir: worldbank

- URB:
Razlaga: Urbano prebivalstvo (% celotne populacije).
Napoved: Za delež urbanega prebivalstva predvidevam, da bodo višje stopnje urbanizacije vodile v višji delež obnovljive energije, ker je urbano prebivalstvo bolj ozaveščeno o zeleni energiji. Prav tako so potrebe po energiji urbanega prebivalstva višje, kar vodi v različne vire pridobivanja energije.
Vir: worldbank

- CO2:
Razlaga: Količina ogljikovega dioksida na prebivalca (v tonah).
Napoved: Večja količina CO2 na prebivalca je povezana z večjim onesnaževanjem in večjo uporabo fosilnih goriv. Pričakujem, da bodo višje ravni CO2 vodile v nižji delež obnovljive energije. Vir: https://ourworldindata.org/grapher/co-emissions-per-capita?tab=table; Global Carbon Budget (2023)

- ELE:
Razlaga: Cena elektrike za ne-gospodinjstva (kilovatne ure; brez davkov in trošarin).
Napoved: Za ceno elektrike za negospodinjstva pričakujem, da bodo višje cene pomenile nižji delež obnovljive energije, saj so zaenkrat cene energije iz obnovljivih virov dražje od neobnovljivih. Hkrati pa višji delež obnovljive energije lahko pomeni nižjo ceno energije iz obnovljivih virov.
Vir: Eurostat; https://ec.europa.eu/eurostat/databrowser/view/nrg_pc_205/default/table?lang=en

- ENG:
Razlaga: Poraba energije na prebivalca (kilovatne ure).
Napoved: Poraba energije na prebivalca je lahko odraz stopnje industrializacije in energetske učinkovitosti države. Pričakujem, da bo višji delež obnovljive energije povezan z manjšo porabo energije. Menim, da višji delež obnovljive energije vodi v višjo ozaveščenost in boljšo energetsko učinkovitost.
Vir: U.S. Energy Information Administration (EIA); https://www.eia.gov/opendata/index.php#bulk-downloads

- GZD:
Razlaga: Gozdno območje (% celotne površine).
Napoved: Pričakujem, da bo višji delež gozdnega območja vodil v višji delež obnovljive energije, saj je biomasa lahko vir obnovljive energije.
Vir: worldbank

library(plm)
podatki_panel <- pdata.frame(podatki, index = c("DRZAVA", "LETO"), drop.inex = FALSE)
head(podatki_panel, 5)
##               DRZAVA LETO   REN      BDP   URB  CO2  ELE      ENG   GZD JOIN
## Austria-2013 Austria 2013 32.67 44299.38 57.34 7.99 0.09 48493.57 46.95 1995
## Austria-2014 Austria 2014 33.55 44245.17 57.53 7.51 0.08 46196.47 46.99 1995
## Austria-2015 Austria 2015 33.50 44195.82 57.72 7.68 0.07 45686.79 47.03 1995
## Austria-2016 Austria 2016 33.37 44590.25 57.91 7.69 0.07 46698.79 47.08 1995
## Austria-2017 Austria 2017 33.14 45281.72 58.09 7.91 0.07 47423.05 47.12 1995
##              drop.inex
## Austria-2013     FALSE
## Austria-2014     FALSE
## Austria-2015     FALSE
## Austria-2016     FALSE
## Austria-2017     FALSE

^ Prikaz prvih 5 vrstic podatkov. Vidmo ustvarjen index države in leta.

podatki_panel_selected <- podatki_panel[, c("REN", "BDP", "URB", "CO2", "ELE", "ENG", "GZD")]
summary(podatki_panel_selected)
##       REN             BDP              URB             CO2        
##  Min.   : 3.49   Min.   :  6701   Min.   :53.33   Min.   : 2.900  
##  1st Qu.:13.86   1st Qu.: 16598   1st Qu.:63.23   1st Qu.: 4.750  
##  Median :18.73   Median : 24785   Median :71.50   Median : 6.135  
##  Mean   :21.93   Mean   : 31859   Mean   :73.30   Mean   : 6.732  
##  3rd Qu.:29.29   3rd Qu.: 43357   3rd Qu.:85.22   3rd Qu.: 8.068  
##  Max.   :66.00   Max.   :110426   Max.   :98.15   Max.   :19.040  
##       ELE               ENG             GZD       
##  Min.   :0.06000   Min.   :18270   Min.   : 1.09  
##  1st Qu.:0.07000   1st Qu.:28185   1st Qu.:22.77  
##  Median :0.08000   Median :35054   Median :34.45  
##  Mean   :0.09815   Mean   :39981   Mean   :34.78  
##  3rd Qu.:0.10750   3rd Qu.:48091   3rd Qu.:39.97  
##  Max.   :0.34000   Max.   :84137   Max.   :73.74

^ Opisna statistika:
- REN: Podatki kažejo, da je povprečni delež energije iz obnovljivih virov 21.93 %. Večina držav ima delež obnovljive energije med 13.86 % in 29.29 %. Minimalna vrednost kaže, da nekatere države porabijo le 3.49 % energije iz obnovljivih virov, medtem ko največja vrednost pokaže na države, ki dosegajo 66.00 %.
- BDP: Povprečni BDP na prebivalca je približno 31.859 €. Mediana je nekoliko nižja, 24.785 €. Razpon vrednosti BDP je precej velik, saj se giblje od 6.701 € do 110.426 €.
- URB: Povprečni delež urbanega prebivalstva je približno 73.30 %. Večina držav ima delež urbanega prebivalstva med 63.23 % in 85.22 %. Minimalna vrednost kaže, da imajo nekatere države le 53.33 % urbanega prebivalstva, medtem ko imajo nekatere 98.15 % urbanega prebivalstva.
- CO2: Povprečna količina CO2 na prebivalca je približno 6.732 ton. Mediana je nekoliko nižja, 6.135 ton. Razpon vrednosti količine CO2 na prebivalca je precej velik, saj se giblje od 2.900 ton do 19.040 ton.
- ELE: Povprečna cena elektrike za negospodinjstva je približno 0.09815 €. Večina držav ima ceno elektrike med 0.07000 € in 0.10750 €.
- ENG: Povprečna poraba energije na prebivalca je približno 39.981 kilovatnih ur. Mediana je nekoliko nižja, 35.054 kilovatnih ur. Razpon vrednosti porabe energije na prebivalca je precej velik, saj se giblje od 18.270 kilovatnih ur do 84.137 kilovatnih ur.
- GZD: Povprečni delež gozdnega območja je približno 34.78 %. Večina držav ima delež gozdnega območja med 22.77 % in 39.97 %. Minimalna vrednost kaže, da imajo nekatere države le 1.09 % gozdnega območja, medtem ko imajo nekatere 73.74 % gozdnega območja.

VIZUALIZACIJA:

library(gplots)
library(plotly)
library(ggplot2)
library(GGally)
library(dplyr)
library(rworldmap)
library(rworldxtra)
map_data <- joinCountryData2Map(podatki_panel, joinCode = "NAME", nameJoinColumn = "DRZAVA")
## 260 codes from your data successfully matched countries in the map
## 10 codes from your data failed to match with a country code in the map
## 217 codes from the map weren't represented in your data
world_map <- getMap()
europe_map <- subset(world_map, REGION == "Europe")


mapCountryData(map_data, nameColumnToPlot="REN", mapTitle="Obnovljiva energija v Evropi", 
               catMethod="fixed", numCats=5, missingCountryCol="gray", addLegend=TRUE, 
               mapRegion="Europe", colourPalette=c("red", "orange", "yellow", "green", "darkgreen"))

^ Vizualni prikaz porazdelitve uporabe obnovljive energije v evropskih državah. Vidimo, da ima Švedska največji povprečni delež tekom desetih let. Finska, Latvija in Avstrija imajo prav tako visok delež. Najnižji delež imajo Irska, Belgija, Nizozemska, Poljska, Lichenstain, Slovaška in Ciper. Upoštevati moramo, da se barva lestvica prične pri 3.49% povprečnega deleža obnovljive energije in ne pri 0. Podatka za Češko ne morem prikazati.

library(dplyr)
library(ggplot2)
library(viridis)

data(podatki, package = "plm")

podatki %>%
  ggplot(aes(x = DRZAVA, y = REN, color = as.factor(LETO))) +
  geom_point() +
  scale_color_viridis_d(name = "LETO") +  # Dodana legenda z imenom "LETO"
  labs(x = "DRZAVA", y = "REN") +
  ggtitle("Gibanje vrednost REN čez države in čas") +
  theme(axis.text.x = element_text(angle = 90))

^ Vizualni prikaz gibanja vrednosti REN za posamezno državo in čez čas. Posamezna pika predstavlja opazovanje deleža energije iz obnovljivih virov za določeno državo za eno leto. Skupaj je 10 opazovanj za posamezno državo. Večja razpršenost točk pomeni večjo variabilnost vrednosti REN tekom desetih let. Vidimo, da imajo nekatere države ,kot so Danska, Estonija, in Švedska visoko razpršenost točk, kar pomeni, da so v deset letnem obdobju močno spreminjale delež obnovljive energije. Nekatere države, kot so Avstrija, Češka, Italija in Romunija, pa imajo nizko razpršenost, kar nam pove, da v desetih letih niso veliko spreminjale deleža obnovljive energije. Pri nekaterih državah lahko celo vidimo, da se je v 10 letih delež obnovljive energije zmanjšal. V primeru Madžarke, je bil delež višji leta 2013 kot 2022.

plotmeans(REN ~ LETO, 
          main="Povprečna vrednost REN čez čas", 
          data= podatki)

^ Prikaz gibanja povprečne vrednosti REN skozi 10 letno obdobje. Vidimo lahko enakomeren trend naraščanja. Vrednost REN je v 10 letih narasla iz cca. 19% na 26%, kar je 7% rast.

my_theme <- theme(axis.text.x = element_text(angle = 90))

podatki %>%
  group_by(LETO) %>%
  summarise(REN_mean = mean(REN)) %>%
  left_join(podatki) %>%
  ggplot(data = ., 
         aes(x = LETO, y = REN)) +
  geom_point() +
  geom_line(aes(x = LETO, y = REN_mean), col = "blue") +
  scale_x_continuous(labels = as.character(podatki$LETO), 
                     breaks = podatki$LETO) +
  labs(x = "LETO", y = "REN") +
  ggtitle("Povprečna vrednost REN čez čas") +
  my_theme

^ Vizulani prikaz gibanja povprečne vrednosti REN skozi 10 letno obdobje. Pike prikazuje dejanske vrednosti REN posamezne države za izbrano leto. Modra črta prikazuje povprečne vrednosti za posamezna leta. Lahko opazimo, da nekaj držav izstopa na višjem delu vrednosti. Vidimo lahko, da je viden dvig povprečne vrednosti med letoma 2019 in 2020. Vidimo pozitiven in enakomeren trend naraščanja povprečnega REN.

ggplot(podatki, aes(x = DRZAVA, y = REN)) +
  geom_boxplot() +
  labs(title = "Porazdelitev kazalnika REN med državami",
       x = "Država",
       y = "Kazalnik REN") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

^ Prikaz porazdelitve kazalnika REN z uporabo boxplot. Boxplot prikaže mediano vrednosti, interkvartilni razpon (razpon med 25. in 75. centilom), iztopajoče vrednosti in razpone vrednosti med max in min. Vidimo, da imajo nekatere države (Avstrija, Romunija) zelo nizko variabilnost vrednost. Nekatere države (Danska, Švedska) pa visoko variabilnost. Črne pike predstavljajo iztopajoče vrednosti in vidimo, da jih ima Estonija kar tri. Vidimo, da ima večina držav vrednost mediane blizu prvega kvartala ali celo minimalne vrednosti (Bulgarija, Češka, Nizozemska).

plot_ly(data = podatki_panel, x = ~LETO, y = ~REN, color = ~DRZAVA, type = "scatter", mode = "lines", hoverinfo = "text",
        text = ~paste("Država: ", DRZAVA, "<br>Leto: ", LETO, "<br>Kazalnik REN: ", REN)) %>%
  layout(title = "Interaktivna vizualizacija trendov v uporabi obnovljivih virov energije",
         xaxis = list(title = "Leto"),
         yaxis = list(title = "Kazalnik REN"),
         hovermode = "closest",
         showlegend = TRUE) %>%
  animation_opts(frame = 100, redraw = TRUE)

^ Prikaz trendov REN po državah skozi leta. Za vsako državo posebaj vidimo, da je količina REN skozi čas rasla. Nekatere države imajo v posameznih letih poskok vrednosti, nekatere države (Danska, Estonija) pa opazne upade. večji premiki, neskladni z trendu, lahko pripišemo različnim dejavnikom, kot so državne subvencije, regulatorne spremembe, naravni dejavniki, nove tehnologije ali politične spremembe. Švedeska ima opazno največjo vrednost. Graf je interaktivne narave. omogoča izbiro časovnega obdobja, posamezne države in izključevanje posameznih držav.

ggpairs(data = podatki_panel[, c("REN", "BDP", "URB", "CO2", "ELE", "ENG", "GZD")])

^ Prikaz odnosov med spremenljivkami. Na diagonali so prikazani histogrami porazdelitev vrednosti spremenljivk. Pod diagonalo so prikazani grafi razpršitev vrednosti posamezne spremenljivke. Zgornji del prikazuje Pearsonove korelacijske koeficiente, ki kvantitativno označujejo moč in smer povezav med spremenljivkami. Vidimo lahko, da so nekatere korelacijo statistično pomembne.
Korelacija med:
CO2 ~ REN; koeficient -0,35, zmerna negativna korelacija, z višanjem deleža REN se manjša količina CO2
GZD ~ REN; koeficient 0,715, visoka pozitivna korelacija, z višanjem deleža GZD, se viša delež REN
BDP ~ URB; koefeicient 0,486, zmerna pozitivna korelacija, države z višjim BDP imajo višji delež urbanizacije
BDP ~ CO2; koeficient 0,536, zmerna pozitivna korelacija, države z višjim BDP imajo višje emisije CO2 na prebivalca
BDP ~ ENG; koeficient 0,618, zmerna pozitivna korelacija, države z višjim BDP porabijo več energije na prebivalca
URB ~ CO2; koeficient 0,220, zmerna pozitivna korelacija, države z višjo stopnjo urbanizacije imajo višje emisije CO2 na prebivalca
URB ~ ENG; koeficient 0,689, visoka pozitivna korelacija, države z višjo stopnjo urbanizacije imajo večjo porabo energije na prebivalca
URB ~ GZD; koeficient -0,221, zmerna negativna korelacija, države z višjo stopnjo urbanizacijo imajo nižji delež gozdnih površin in obratno
ENG ~ CO2; koeficient 0,543, zmerna pozitivna korelacija, višja poraba energije vodi v več emisij CO2 na prebivalca
GZD ~ ELE; koeficient -0,260, zmerna negativna korelacija, države ki imajo išjo ceno elektrike imajo manjši delež gozdnih površin

my_colors <- c("#1F77B4", "#FF7F0E", "#2CA02C", "#D62728", "#9467BD", "#8C564B", "#E377C2", "#7F7F7F",
               "#BCBD22", "#17BECF", "#1A1A1A", "#7FDBFF", "#2ECC40", "#FFDC00", "#FF851B", "#B10DC9",
               "#F012BE", "#FF4136", "#85144B", "#AAAAAA", "#DDDDDD", "#333333", "#0074D9", "#3D9970",
               "#FF4136", "#B10DC9", "#FF851B")


ggplot(data = podatki, aes(x = BDP, y = REN)) +
  geom_point(aes(color = factor(DRZAVA))) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = my_colors) +  # Uporaba nove barvne palete
  labs(x = "BDP",
       y = "REN",
       color = "DRZAVA")

^ Grafični prikaz razmerja med BDP in REN. Vsaka točka predstavlja vrednost leta posamezne države. Barva točke je določena glede na državo. regresijska črta je rahlo padajoča. Naraščanje BDP ne vodi v sorazmerno naraščanje REN. Vidimo, da nekatere države izstopajo. Irska in Luksemburg z visokim BDP vednar nisko stopnjo REN. Ter Finska in Švedska z visoko REN in nizkim BDP.

ggplot(data = podatki,
       aes(x = URB, y = REN)) +
  geom_point(aes(color = factor(DRZAVA))) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = my_colors) +
  labs(x = "URB",
       y = "REN",
       color = "DRZAVA")

^ Prikaz razmerja med URB in REN. Regresijska črta je rahlo padajoča. Naraščanje URB ne vodi v sorazmerno naraščanje BDP. Med spremenljivkami je rahla negativna liearna povezava. Vidimo, da so nekatere točke držav navpične, kar pomeni, da se z višanjem obnovljive energije ni spreminjala stopnja urbanizacije.

ggplot(data = podatki,
       aes(x = CO2, y = REN)) +
  geom_point(aes(color = factor(DRZAVA))) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = my_colors) +
  labs(x = "CO2",
       y = "REN",
       color = "DRZAVA")

^ Prikaz razmerja med CO2 in REN. Regresijska črta je močno padajoča. Med spremenljivkama je izrazita negativna linearna povezava. Naraščanje CO2 na prebivalaca vodi v izrazito nižanje stopnje obnovljive energije. Vidimo, da Luksemburg izstopa z visoko količino CO2 in nizko stopnjo obnovljive energije.

ggplot(data = podatki,
       aes(x = ELE, y = REN)) +
  geom_point(aes(color = factor(DRZAVA))) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = my_colors) +
  labs(x = "ELE",
       y = "REN",
       color = "DRZAVA")

^ Prikaz razmerja med ELE in REN. Regresijska črta je padajoča. Med spremenljivkama je izrazita negativna linearna povezava. Naraščanje EL ne vodi v sorazmerno naraščanje BDP. Vidimo, da je veliko enot navpičnih, kar pomeni, da se z cena elektrike skozi leta ni drastičnos spreminjala, medtem ko se je spreminjal delež obnovljive energije.

ggplot(data = podatki,
       aes(x = ENG, y = REN)) +
  geom_point(aes(color = factor(DRZAVA))) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = my_colors) +
  labs(x = "ENG",
       y = "REN",
       color = "DRZAVA")

^ Prikaz razmerja med URB in REN. Regresijska črta je rahlo padajoča. Med spremenljivkama je negativna linearna povezava. Naraščanje ENG ne vodi v sorazmerno naraščanje BDP. Vidimo, da imajo države z višjo porabo energije na prebivalca nižje stopnje obnovljive energije.

ggplot(data = podatki,
       aes(x = GZD, y = REN)) +
  geom_point(aes(color = factor(DRZAVA))) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = my_colors) +
  labs(x = "GZD",
       y = "REN",
       color = "DRZAVA")

^ Prikaz razmerja med GZD in REN. Regresijska črta je močno naraščujoča. Med spremenljivkama je izrazita pozitivna linearna povezava. Države z višjim deležem gozda imajo večji delež obnovljive energije. Vidimo, da se delež gozda pri večini držav skozi leta ni spreminjal.

REGRESIJA:

library(lmtest)
library(car)
library(plm)
panel_model <- plm(REN ~ BDP + URB + CO2 + ELE + ENG + GZD, 
                   data = podatki_panel, 
                   index = c("DRZAVA", "LETO"), model = "random")

^ Prvi panelni model. Za prvi model sem izbral naključne učnike, rezultati testov se ne spremenijo glede na uporabo FE in RE. Panel_model bom uporabil za preverjanje predpostavk panelne regresije.

plmtest(panel_model, effect = "individual")
## 
##  Lagrange Multiplier Test - (Honda)
## 
## data:  REN ~ BDP + URB + CO2 + ELE + ENG + GZD
## normal = 30.5, p-value < 2.2e-16
## alternative hypothesis: significant effects

^ Test za preverjanje heteroskedastičnostiv napakah modela. Testna statistika je statistično pomembna (p<0,05), kar nakazuje na prisotnost heteroskedastičnosti. Ker v modelu ne želimo heteroskedastičnosti, bom v nadaljevanju prilagodil panelni model z uporabo robustnih standardnih napak.

bptest(panel_model)
## 
##  studentized Breusch-Pagan test
## 
## data:  panel_model
## BP = 16.353, df = 6, p-value = 0.01198

^ Dodaten test za preverjanje heteroskedastičnosti. Test ponovno kaže na statistično pomembno heteroskedastičnost v panelnem modelu (p<0,05).

shapiro.test(residuals(panel_model))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(panel_model)
## W = 0.97601, p-value = 0.0001639

^ Test normalnosti pokaže, da so ostanki porazdeljeni blizu normalne porazdelitve, saj je vrednost testne statistike W = 0.97601. Kljub temu pa je p-vrednost manjša od 0.05, kar kaže, da se ostanki ne porazdeljujejo normalno. Zato bom nadaljnjo analizo prilagodil z uporabo robustnih standardnih napak.

vif(panel_model)
##      BDP      URB      CO2      ELE      ENG      GZD 
## 1.282652 1.376616 1.848615 1.121416 1.920317 1.045510

^ Test za preverjanje prisotnosti multikolinearnosti med spremenljivkami kaže, da so vse vrednosti variacijskega faktorja (VIF) pod 5, kar nakazuje na nizko stopnjo multikolinearnosti. To pomeni, da pojasnjevalne spremenljivke v modelu niso močno povezane med seboj, kar omogoča zanesljivo ocenjevanje parametrov modela.

FE_panel_model <- plm(REN ~ BDP + URB + CO2 + ELE + ENG + GZD, data = podatki,  
                      index = c("DRZAVA", "LETO"), 
                      effect = "individual", model = "within")

^ Panelni model z uporabo fiksnih učinkov. Fiksni učinki v panelnih modelih obravnavajo stalne značilnosti posameznih enot. Pri uporabi fiksnih učinkov se ocenjuje vpliv spremenljivk na odvisno spremenljivko znotraj vsake enote posebej, kar omogoča upoštevanje stalnih značilnosti enot, ki se ne spreminjajo s časom.

RE_panel_model <- plm(REN ~ BDP + URB + CO2 + ELE + ENG + GZD, data = podatki,  
                      index = c("DRZAVA", "LETO"), 
                      effect = "individual", model = "random")

^ Panelni model z uporabo naključnimh učinkov. Naključni učinki v panelnih modelih obravnavajo nespremenljive značilnosti enot kot vzorce iz naključne populacije. Pri uporabi naključnih učinkov se ocenjuje vpliv spremenljivk na odvisno spremenljivko z upoštevanjem naključnih variacij v učinkih med enotami. To omogoča upoštevanje neopazovanih heterogenosti med enotami, ki se lahko spreminjajo s časom.

POOLED_panel_model <- plm(REN ~ BDP + URB + CO2 + ELE + ENG + GZD, data = podatki, 
                      index = c("DRZAVA", "LETO"), 
                      effect = "individual", model = "pooling")

^ Panelni model z metodo združevanja. Metoda združevanja v analizi panelnih podatkov predpostavlja, da so učinki pojasnjevalnih spremenljivk enaki za vse enote v panelu, kar omogoča ocenjevanje povprečnih učinkov preko vseh enot brez upoštevanja individualnih razlik med njimi. Ta pristop se običajno uporablja, kadar ni pričakovati pomembnih individualnih razlik med enotami ali ko so te razlike nepomembne za namen analize.

phtest(FE_panel_model, RE_panel_model)
## 
##  Hausman Test
## 
## data:  REN ~ BDP + URB + CO2 + ELE + ENG + GZD
## chisq = 26.229, df = 6, p-value = 0.0002018
## alternative hypothesis: one model is inconsistent

^ Test za primerjo razlik med ocenami parametrov med modeloma FE in RE. Ker je p-vrednost manjša od 0,05, to kaže na to, da je eden izmed modelov neskladen z podatki. Glede na rezultate testa je bolje izbrati pristop modela z fiksnimi učinki (FE).

plmtest(POOLED_panel_model, effect = "individual", type = c("bp"))
## 
##  Lagrange Multiplier Test - (Breusch-Pagan)
## 
## data:  REN ~ BDP + URB + CO2 + ELE + ENG + GZD
## chisq = 930.25, df = 1, p-value < 2.2e-16
## alternative hypothesis: significant effects

^ Test Lagrange Multiplier Breusch-Pagan pomaga pri izbiri med modelom naključnih učinkov in preprostim regresijskim modelom OLS. Rezultati testa kažejo, da obstajajo pomembne razlike med državami (p < 0,05). Zato združena regresija OLS ni ustrezna, boljša izbira je model naključnih učinkov.

FE_panel_model_RBS <- plm(REN ~ BDP + URB + CO2 + ELE + ENG + GZD, 
                      data = podatki,  
                      index = c("DRZAVA", "LETO"), 
                      effect = "individual", 
                      model = "within")
robust_se <- coeftest(FE_panel_model_RBS, vcov = vcovHC(FE_panel_model_RBS, type = "HC1"))
print(robust_se)
## 
## t test of coefficients:
## 
##        Estimate  Std. Error t value  Pr(>|t|)    
## BDP  1.5817e-04  1.1376e-04  1.3903 0.1657372    
## URB  1.1750e+00  3.7645e-01  3.1214 0.0020234 ** 
## CO2 -1.3608e+00  3.5903e-01 -3.7903 0.0001909 ***
## ELE  3.9733e+00  3.6430e+00  1.0907 0.2765159    
## ENG -8.8117e-06  5.8176e-05 -0.1515 0.8797369    
## GZD  3.1002e-01  6.3211e-01  0.4905 0.6242615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(FE_panel_model_RBS)
## 
## Model Formula: REN ~ BDP + URB + CO2 + ELE + ENG + GZD
## 
## Coefficients:
##         BDP         URB         CO2         ELE         ENG         GZD 
##  1.5817e-04  1.1750e+00 -1.3608e+00  3.9733e+00 -8.8117e-06  3.1002e-01

^ Panelna regresija z uporabo fiksnih efektov in uporabo robustnih standardih napak PRAVILNI MODEL
Od 6 pojasnjevalnih spremenljivk, sta statistično značilni samo URB in CO2. Spremenljivke BDP, ELE, ENG in GZD nimajo statistično pomembnega vpliva na odvisno spremenljivko delež obnovljive energije.

- URB: Če se delež urbanizacije poveča za 1 odstotno točko, se delež obnovljive energije poveča za 1.17 odstotne točke (p=0,002), ob predpostavki, da so vse ostale pojasnjevalne spremenljivke nespremenjene.
- CO2: Če se količina ogljikovega dioksida na prebivalca poveča za eno tono, se delež obnovljive energije zniža za 1,36 odstotne točke (p=0,0002), ob predpostavki, da so vse ostale pojasnjevalne spremenljivke nespremenjene.

summary(FE_panel_model)
## Oneway (individual) effect Within Model
## 
## Call:
## plm(formula = REN ~ BDP + URB + CO2 + ELE + ENG + GZD, data = podatki, 
##     effect = "individual", model = "within", index = c("DRZAVA", 
##         "LETO"))
## 
## Balanced Panel: n = 27, T = 10, N = 270
## 
## Residuals:
##      Min.   1st Qu.    Median   3rd Qu.      Max. 
## -7.550312 -0.996862 -0.094606  0.977465  7.260825 
## 
## Coefficients:
##        Estimate  Std. Error t-value  Pr(>|t|)    
## BDP  1.5817e-04  3.9903e-05  3.9638 9.767e-05 ***
## URB  1.1750e+00  2.1846e-01  5.3788 1.797e-07 ***
## CO2 -1.3608e+00  1.9889e-01 -6.8423 6.597e-11 ***
## ELE  3.9733e+00  3.5669e+00  1.1140    0.2664    
## ENG -8.8117e-06  5.8337e-05 -0.1510    0.8801    
## GZD  3.1002e-01  5.8479e-01  0.5301    0.5965    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Total Sum of Squares:    1958.8
## Residual Sum of Squares: 855.5
## R-Squared:      0.56325
## Adj. R-Squared: 0.50429
## F-statistic: 50.9418 on 6 and 237 DF, p-value: < 2.22e-16

^ Panelna regresija z uporabo fiksnih efektov. ZA PRIMERJAVO
V primerjavi z panelno regresijo z uporabo fiksnih učinkov in robustnimi standardnimi napakami, je postala spremenljivka BDP statistično značilna in spremenljivka URB je pridobila na statistični moči.

- BDP: Če se bruto domači proizvod na prebivalca poveča za 1 Eur, se delež obnovljive energije poveča za 0,00016 odstotne točke (p=0,0,0001), ob predpostavki, da so vse ostale pojasnjevalne spremenljivke nespremenjene.
- URB: Če se delež urbanizacije poveča za 1 odstotno točko, se delež obnovljive energije poveča za 1.17 odstotne točke (p<0,0001), ob predpostavki, da so vse ostale pojasnjevalne spremenljivke nespremenjene.
- CO2: Če se količina ogljivkovega dioksida na prebivalca poveča za eno tono, se delež obnovljive energije zniža za 1,36 odstotne točke (p<0,0001), ob predpostavki, da so vse ostale pojasnjevalne spremenljivke nespremenjene.

RE_panel_model_RBS <- plm(REN ~ BDP + URB + CO2 + ELE + ENG + GZD, 
                      data = podatki,  
                      index = c("DRZAVA", "LETO"), 
                      effect = "individual", 
                      model = "random")
robust_se <- coeftest(RE_panel_model_RBS, vcov = vcovHC(RE_panel_model_RBS, type = "HC1"))
print(robust_se)
## 
## t test of coefficients:
## 
##                Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept) -1.1152e+01  1.0833e+01 -1.0295 0.3042010    
## BDP          1.6733e-04  8.0446e-05  2.0801 0.0384881 *  
## URB          2.6367e-01  1.4432e-01  1.8270 0.0688367 .  
## CO2         -1.5773e+00  4.5024e-01 -3.5033 0.0005399 ***
## ELE          8.9246e+00  3.7109e+00  2.4050 0.0168636 *  
## ENG         -6.0107e-05  6.6790e-05 -0.8999 0.3689807    
## GZD          5.9151e-01  6.8634e-02  8.6182 6.441e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print(RE_panel_model_RBS)
## 
## Model Formula: REN ~ BDP + URB + CO2 + ELE + ENG + GZD
## 
## Coefficients:
## (Intercept)         BDP         URB         CO2         ELE         ENG 
## -1.1152e+01  1.6733e-04  2.6367e-01 -1.5773e+00  8.9246e+00 -6.0107e-05 
##         GZD 
##  5.9151e-01

^ izračun panelne regresije z naključnimi učinki in uporabo robustnih standardnih napak. ZA PRIMERJAVO

TRENDI:

arima_model <- arima(podatki$REN, order = c(1, 0, 0))
future_values <- predict(arima_model, n.ahead = 10)
print(future_values)
## $pred
## Time Series:
## Start = 271 
## End = 280 
## Frequency = 1 
##  [1] 63.34637 60.85815 58.52503 56.33733 54.28599 52.36252 50.55894 48.86778
##  [9] 47.28204 45.79513
## 
## $se
## Time Series:
## Start = 271 
## End = 280 
## Frequency = 1 
##  [1]  4.768731  6.537202  7.766231  8.704625  9.453035 10.065186 10.574166
##  [8] 11.002237 11.365295 11.675180

^ Model napovedovanja bodočih vrednosti z uporabo ARIMA modela. Test napoveduje vrednosti za naslednjih 10 let. Napoved kaže na pričakovano zmanjšanje deleža obnovljive energije v prihodnjih 10 letih. upad v desetih letih naj bi bil iz 63,35% na 45,79%. Standardne napake so visoke, kar kaže na razmeroma visoko negotovost v napovedi.

library(forecast)
year_start <- 2013
month_start <- 1
year_end <- 2022
month_end <- 12

time_series <- ts(podatki$REN, start = c(year_start, month_start), end = c(year_end, month_end), frequency = 12)
holt_winters_model <- HoltWinters(time_series)
forecast_values <- forecast(holt_winters_model, h = 100, level = c(95))
plot(forecast_values)

^ Model za napovedovanje prihodnjih vrednosti REN. Projekcija napoveduje naraščanje vrednosti deleža obnovljive energije v prihodnosti. to nakazuje na pozitiven trend rasti, ki je v navzkrižju z ARIMA modelom, ki pa napoveduje upadanje.

ZAKLJUČEK: