Wstęp do projektu zaliczeniowego z przedmiotu “Prognozowanie i Symulacje”

W ramach niniejszego projektu analizie poddane zostało zjawisko gospodarcze dotyczące zależności pomiędzy nakładami inwestycyjnymi (Y) a eksportem, PKB i inflacją (X) w Polsce. Badany okres obejmuje lata 2010-2020, a prognoza właściwa została przygotowana na rok 2021.

Wybór zmiennych niezależnych, takich jak eksport, PKB i inflacja, jest uzasadniony ich istotnym wpływem na nakłady inwestycyjne. Wzrost eksportu często prowadzi do zwiększenia inwestycji w infrastrukturę produkcyjną. Analogicznie, zmiany w PKB wpływają na zdolność firm do inwestowania, a inflacja ma bezpośredni wpływ na koszty inwestycji.

Analiza ekonomiczna Polski w latach 2010-2020 wskazuje na dynamiczny rozwój gospodarki. Realne PKB Polski rosło średnio o 3,2% rocznie, osiągając w 2021 roku wartość 677,5 miliarda dolarów. Wzrost gospodarczy był wspierany przez stabilny eksport i umiarkowaną inflację, choć w ostatnich latach analizowanego okresu wzrost inflacji oraz wpływ globalnych wydarzeń, takich jak: pandemia COVID-19, przerwane łańcuchy dostaw i recesja wśród czołowej gospodarki Europy i głównego partnera Polski - Niemiec stanowiły zauważalne wyzwanie. W literaturze przedmiotu, w tym raportach i analizach dostępnych w bazach danych OECD i World Bank, dostarczyły wartościowych informacji do przygotowania prognoz i analiz trendów gospodarczych w Polsce.

Celem ów projektu jest zbudowanie modelu przyczynowo-skutkowego, który pozwoli na lepsze zrozumienie mechanizmów wpływających na decyzje inwestycyjne oraz przygotowanie prognozy na rok 2021, bazując na zgromadzonych danych historycznych oraz analizie ekonomicznej.

Projekt wykonany został w formacie HMTL, a jego układ został przygotowany tak, aby prowadzić czytelnika za przysłowiową rekę. Kolejne tytuły zawierają wstęp do każdego z rozdziałów, zaś opisy, do wykresów oraz modeli zostały umieszczone poniżej tychże elementów.

Wyświetlenie danych i informacji o nich

Budowa modelu przyczynowo-skutkowego rozpoczęła się od wyszukania danych i wprowadzenia ich do jednolitego formatu, a kolejno do pliku, który wczytany został do programu R.

dane <- read_excel('NOWE.xlsx')

Porządkowanie danych

Pomimo wcześniejszego zebrania i przedstawienia danych w jednym pliku konieczne było ich przetworzenie i przygotowanie do dalszej eksploracji oraz prowadzenia badań.

dane <- dane %>%
  mutate(rok = make_date(year = rok, month = 1, day = 1))

ustawienie_kolumn <- c('rok', 'inf', 'Eks', 'PKB', 'ni')
dane <- dane %>%
  select(all_of(ustawienie_kolumn))

Wizualizacje rozkładów zmiennych oraz ich zachowania w czasie

Każda ze zmiennych przedstawiona została na dwa sposoby. Po pierwsze zobrazowany został rozkład każdej z nich za pomocą wykresu typ histogram. Drugim rodzajem był prosty wykres liniowy, którego celem było przedstawienie charakteru zmiennej w czasie.

ggplot(dane, aes(x = PKB)) +
  geom_histogram(bins = 5, fill = "blue", color = "white") + 
  theme_minimal() + 
  labs(x = 'PKB', y = 'Rozkład', title = 'Rozkład danych o PKB')

ggplot(dane, aes(x = rok, y = PKB)) +
  geom_line(group = 1, color = "red") + 
  theme_minimal() + # Styl wykresu
  labs(x = 'Rok', y = 'PKB', title = 'Zmiana wartości PKB w czasie')

W przypadku danych o Produkcie Krajowym Brutto (PKB) histogram wskazuje, że najwięcej obserwacji z lat 2010-2020 przyjmowało wartość zbliżoną do 1,5 biliona złotych. Rozkład wartości przyjął ksztłat prawoskośnego równego, gdzie znacząco różniła się tylko jedna, wspomniana już grupa wartości. Linia trendu na wykresie czasowym pokazuje systematyczny wzrost wartości PKB.

ggplot(dane, aes(x = Eks)) +
  geom_histogram(bins = 5, fill = "green", color = "white") + 
  theme_minimal() + 
  labs(x = 'Eks', y = 'Rozkład', title = 'Rozkład danych o Eksporcie')

ggplot(dane, aes(x = rok, y = PKB)) +
  geom_line(group = 1, color = "purple") + 
  theme_minimal() + # Styl wykresu
  labs(x = 'Rok', y = 'Eks', title = 'Zmiana wartości Eksportu w czasie')

Druga para wykresów odnosi się do zmiennej ‘Eksport’. Rozkład tej zmiennej również jest prawoskośny, jednakże w tym przypadku dane są w większyms stopniu zróżnicowane niż dla wcześniejszej zmiennej. Wykres czasowy pokazuje tak naprawdę prostą linię w okresie 2010-2020. Takie zjawisko uznać należy za pozytywne, gdyż śwsiadczy o rosnącej konkurencyjności polskiej gospodarki na rynkach zagranicznych.

ggplot(dane, aes(x = ni)) +
  geom_histogram(bins = 5, fill = "pink", color = "white") + 
  theme_minimal() + 
  labs(x = 'ni', y = 'Rozkład', title = 'Rozkład danych o nakładach inwestycyjnych')

ggplot(dane, aes(x = rok, y = ni)) +
  geom_line(group = 1, color = "orange") + 
  theme_minimal() + # Styl wykresu
  labs(x = 'Rok', y = 'ni', title = 'Zmiana wartości nakładów inwestycyjnych w czasie')

Przetworzenie danych o nakładach inwestycyjnych pozwoliło wejrzeć w ich specyfikę i stwierdzić, że zmienność danych była bardziej zauważalna niż w przedstawionych powyżej dwóch zmiennych. Zależność tej zmiennej w czasie ujawniła trend wzrostowy w badanym okresie. Co ciekawe, w przypadku tej zmiennej zauważalne są wahania, które mogły wynikać z nieregularnymi zmianami w polityce inwestycyjnej.

ggplot(dane, aes(x = inf)) +
  geom_histogram(bins = 5, fill = "gray", color = "white") + 
  theme_minimal() + 
  labs(x = 'inf', y = 'Rozkład', title = 'Rozkład danych o inflacji')

ggplot(dane, aes(x = rok, y = inf)) +
  geom_line(group = 1, color = "red") + 
  theme_minimal() + # Styl wykresu
  labs(x = 'Rok', y = 'inf', title = 'Zmiana wartości inflacji w czasie')

Ostatnia para wykresów dotyczy zmiennej obajśniającej - inflacji. Rozkład danych sugeruje relatywnie stabilny poziom inflacji z medianą w okolicach 102. Zestawienie i zobrazowanie inflacji względem czasu, pozwala stwierdzić, że ta zmienna w zwyczaju ma dążenie do charakteru nieliniowego. Znaczące wahania w latach 2010-2020 wynikały ze zmian ekonomicznych bezpośrednio w Polsce.

dane <- dane %>%
  mutate(across(where(is.numeric), list(ln = ~log(.)), .names = "ln_{.col}")) %>%  
  select(-where(is.numeric), everything())
dane

## # A tibble: 11 × 9
##    rok          inf      Eks     PKB      ni ln_inf ln_Eks ln_PKB ln_ni
##    <date>     <dbl>    <dbl>   <dbl>   <dbl>  <dbl>  <dbl>  <dbl> <dbl>
##  1 2010-01-01 103.   481058. 1434368 114075.   4.63   13.1   14.2  11.6
##  2 2011-01-01 104.   558739  1553641 130689.   4.65   13.2   14.3  11.8
##  3 2012-01-01 104.   603419. 1612739 128097.   4.64   13.3   14.3  11.8
##  4 2013-01-01 101.   647879. 1630126 132956.   4.61   13.4   14.3  11.8
##  5 2014-01-01 100    693472. 1700552 150294.   4.61   13.4   14.3  11.9
##  6 2015-01-01  99.1  750836. 1798471 167157.   4.60   13.5   14.4  12.0
##  7 2016-01-01  99.4  803478. 1853205 149087.   4.60   13.6   14.4  11.9
##  8 2017-01-01 102    882620. 1982794 159098.   4.62   13.7   14.5  12.0
##  9 2018-01-01 102.   951324. 2126506 181487.   4.62   13.8   14.6  12.1
## 10 2019-01-01 102.  1023591. 2288492 199895.   4.63   13.8   14.6  12.2
## 11 2020-01-01 103.  1062514. 2337672 184295.   4.64   13.9   14.7  12.1

Badanie stacjonarności - test ADF

Problematyka niniejszego projektu sprowadza się do prognozowania szeregów czasowych. Kluczowym krokiem w analizie takiego zjawiska jest badanie stacjonarności. Stacjonarność odnosi się do właściwości procesu, w którym jego statystyki, takie jak średnia, wariancja i autokorelacja, pozostają niezmienne w czasie. Istnieje wiele testów na stacjonarność, najpopularniejszym jest test rozszerzony Augmented Dickeya-Fullera (ADF). W ramach tego projektu wykorzystana została ta metoda statystyczna.

test_adf <- lapply(dane[sapply(dane, is.numeric)], adf.test)

## Warning in FUN(X[[i]], ...): p-value smaller than printed p-value

statystyki_tabela <- data.frame(
  Zmienna = c('Inflacja', 'Eksport', 'PKB', 'Nakłady Inwestycyjne'),
  Statystka_ADF = sapply(test_adf, function(x) x$statistic),
  Wartość_p = sapply(test_adf, function(x) x$p.value)
)

Dla zmiennej Eksport wartość p wyniosła 0.2643 i jest wyższa od poziomu istotności 5%, zatem nie ma podstaw do odrzucenia hipotezy zerowej i stwierdzamy, że ta zmienna nie jest stacjonarna.

Dla zmiennej PKB wartość p wyniosła 0.9518 i jest wyższa od poziomu istotności 5%, zatem nie ma podstaw do odrzucenia hipotezy zerowej i stwierdzamy, że ta zmienna nie jest stacjonarna.

Dla zmiennej Nakłady Inwestycyjne wartość p wyniosła 0.7904 i jest wyższa od poziomu istotności 5%, zatem nie ma podstaw do odrzucenia hipotezy zerowej i stwierdzamy, że ta zmienna nie jest stacjonarna.

Dla zmiennej Inflacja wartość p wyniosła 2e-12 i jest niższa od poziomu istotności 5%, zatem odrzucamy hipotezę zerową i stwierdzamy, że ta zmienna jest stacjonarna.

Wykresy ACF i PACF dla zmiennych

Informacje dotyczące stacjonarności zmiennych w opracowywanym modelu stanowią podstawę podejścia do właściwych prognoz, to oznacza dla każdej z determinant nakładów inwestycyjnych osobno. Jedną z wybranych metod prognozowania jest model ARIMA. W celu uzyskania jak najlepszych rezultatów oraz jakości pracy, wykorzystany został model Auto ARIMA. Niemniej jednak dla każdej zmiennej wykonano sprawdzenie funkcji autokorelacji i autokorelacji cząstkowej.

dane_eks <- dane %>% select(rok, ln_Eks)

acf_result <- Acf(dane_eks$ln_Eks, lag.max=5, plot=FALSE)
pacf_result <- Pacf(dane_eks$ln_Eks, lag.max=5, plot=FALSE)

acf_df <- data.frame(lag = acf_result$lag[-1], acf = acf_result$acf[-1]) 
pacf_df <- data.frame(lag = pacf_result$lag[-1], pacf = pacf_result$acf[-1])

p1 <- ggplot(acf_df, aes(x = lag, y = acf)) +
  geom_point() +
  geom_line() +
  labs(title = "Autokorelacja funkcji (ACF)", x = "Opóźnienie", y = "ACF") +
  theme_minimal()

p2 <- ggplot(pacf_df, aes(x = lag, y = pacf)) +
  geom_point() +
  geom_line() +
  labs(title = "Częściowa autokorelacja funkcji (PACF)", x = "Opóźnienie", y = "PACF") +
  theme_minimal()

grid.arrange(p1, p2, nrow=2)

W przypadku zmiennej ‘Eksport’ występuje autokorelacja pierwszego rzędu, co przedstawiają dwa powyższe wykresy. Zjawisko zanika po tym opóźnieniu. Odpowiednim modelem dla zmiennej eksport jest zatem część AR(1) w ARMA.

dane_pkb <- dane %>% select(rok, ln_PKB)

acf_result <- Acf(dane_pkb$ln_PKB, lag.max=5, plot=FALSE)
pacf_result <- Pacf(dane_pkb$ln_PKB, lag.max=5, plot=FALSE)

acf_df <- data.frame(lag = acf_result$lag[-1], acf = acf_result$acf[-1]) 
pacf_df <- data.frame(lag = pacf_result$lag[-1], pacf = pacf_result$acf[-1])

p1 <- ggplot(acf_df, aes(x = lag, y = acf)) +
  geom_point() +
  geom_line() +
  labs(title = "Autokorelacja funkcji (ACF)", x = "Opóźnienie", y = "ACF") +
  theme_minimal()

p2 <- ggplot(pacf_df, aes(x = lag, y = pacf)) +
  geom_point() +
  geom_line() +
  labs(title = "Częściowa autokorelacja funkcji (PACF)", x = "Opóźnienie", y = "PACF") +
  theme_minimal()

grid.arrange(p1, p2, nrow=2)

Również w przypadku zmiennej ‘PKB’ wykresy ACF i PACF wskazują na istotną korelację dla pierwszego opóźnienia. Podobnie, jak wyżej należy spodziewać się, że odpowiednim podejściem modelowania zmiennej ‘PKB’ będzie element AR(1).

dane_ni <- dane %>% select(rok, ln_ni)

acf_result <- Acf(dane_ni$ln_ni, lag.max=5, plot=FALSE)
pacf_result <- Pacf(dane_ni$ln_ni, lag.max=5, plot=FALSE)

acf_df <- data.frame(lag = acf_result$lag[-1], acf = acf_result$acf[-1]) 
pacf_df <- data.frame(lag = pacf_result$lag[-1], pacf = pacf_result$acf[-1])

p1 <- ggplot(acf_df, aes(x = lag, y = acf)) +
  geom_point() +
  geom_line() +
  labs(title = "Autokorelacja funkcji (ACF)", x = "Opóźnienie", y = "ACF") +
  theme_minimal()

p2 <- ggplot(pacf_df, aes(x = lag, y = pacf)) +
  geom_point() +
  geom_line() +
  labs(title = "Częściowa autokorelacja funkcji (PACF)", x = "Opóźnienie", y = "PACF") +
  theme_minimal()

grid.arrange(p1, p2, nrow=2)

W przypadku zmiennej zależnej, to jest nakłady inwestycyjne, istotna okazała się analogicznie dodatnia korelacja pierwszego rzędu. Wykresy sugerują, że odpowiednim podejściem do modelowania będzie przyjęcie elementu AR(1).

dane_inf <- dane %>% select(rok, ln_inf)

acf_result <- Acf(dane_inf$ln_inf, lag.max=5, plot=FALSE)
pacf_result <- Pacf(dane_inf$ln_inf, lag.max=5, plot=FALSE)

acf_df <- data.frame(lag = acf_result$lag[-1], acf = acf_result$acf[-1]) 
pacf_df <- data.frame(lag = pacf_result$lag[-1], pacf = pacf_result$acf[-1])

p1 <- ggplot(acf_df, aes(x = lag, y = acf)) +
  geom_point() +
  geom_line() +
  labs(title = "Autokorelacja funkcji (ACF)", x = "Opóźnienie", y = "ACF") +
  theme_minimal()

p2 <- ggplot(pacf_df, aes(x = lag, y = pacf)) +
  geom_point() +
  geom_line() +
  labs(title = "Częściowa autokorelacja funkcji (PACF)", x = "Opóźnienie", y = "PACF") +
  theme_minimal()

grid.arrange(p1, p2, nrow=2)

Tu w końcu nie wiem o co chodzi :(

ARIMA

Badanie autokorelacji zostało wykonane, a wyniki dotyczące opóźnień zostaną wykorzystane do porównania względem tych zaproponowanych przez model Auto ARIMA.

dane_list <- list(
  dane_eks = "ln_Eks",
  dane_pkb = "ln_PKB",
  dane_ni = "ln_ni",
  dane_inf = "ln_inf"
)

for (dane_name in names(dane_list)) {
  column_name <- dane_list[[dane_name]]
  data <- get(dane_name)  
  
  model <- auto.arima(data[[column_name]])
  print(summary(model))
  checkresiduals(model)
}

## Series: data[[column_name]] 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       0.0792
## s.e.  0.0086
## 
## sigma^2 = 0.0008323:  log likelihood = 21.91
## AIC=-39.82   AICc=-38.1   BIC=-39.21
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.001182227 0.02609547 0.01670845 0.01072918 0.1239864 0.2108576
##                   ACF1
## Training set 0.1142336

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0) with drift
## Q* = 0.93241, df = 3, p-value = 0.8176
## 
## Model df: 0.   Total lags used: 3
## 
## Series: data[[column_name]] 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       0.0488
## s.e.  0.0072
## 
## sigma^2 = 0.0005956:  log likelihood = 23.66
## AIC=-43.31   AICc=-41.6   BIC=-42.71
## 
## Training set error measures:
##                       ME       RMSE       MAE         MPE      MAPE      MASE
## Training set 0.001284308 0.02207483 0.0199428 0.008871638 0.1382242 0.4083032
##                    ACF1
## Training set 0.08143965

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0) with drift
## Q* = 2.4891, df = 3, p-value = 0.4773
## 
## Model df: 0.   Total lags used: 3
## 
## Series: data[[column_name]] 
## ARIMA(0,1,0) 
## 
## sigma^2 = 0.009733:  log likelihood = 8.98
## AIC=-15.96   AICc=-15.46   BIC=-15.66
## 
## Training set error measures:
##                      ME       RMSE        MAE       MPE      MAPE      MASE
## Training set 0.04466563 0.09406265 0.08388293 0.3729739 0.7004169 0.9207103
##                    ACF1
## Training set -0.2323088

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,0)
## Q* = 5.0089, df = 3, p-value = 0.1711
## 
## Model df: 0.   Total lags used: 3
## 
## Series: data[[column_name]] 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.6151  4.6252
## s.e.  0.2243  0.0090
## 
## sigma^2 = 0.0001968:  log likelihood = 32.19
## AIC=-58.38   AICc=-54.95   BIC=-57.19
## 
## Training set error measures:
##                         ME       RMSE       MAE          MPE      MAPE
## Training set -0.0004239144 0.01269062 0.0110754 -0.009958408 0.2397021
##                   MASE      ACF1
## Training set 0.9393142 0.3120976

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 6.8262, df = 3, p-value = 0.07765
## 
## Model df: 1.   Total lags used: 4

Projekt zespołowy - Prognozowanie i Symulacje