Forecasting and Simulations
Łebed, Karpicka, Itrych
General information
Project objective
The primary objective of this project is to model the unemployment rate in Poland over the period 2010-2023 based on selected economic and social indicators: GDP per capita, voivodeship budget revenue and expenditure per capita, university graduates per 10,000 inhabitants, and the inflation rate. The secondary goal of the study is to use the constructed model to forecast unemployment rates in Poland up to the year 2030. All data used in this analysis were obtained from the official database of Statistics Poland (Główny Urząd Statystyczny, GUS).
Explanation of variables
employment= β0 + β1 * GDPpc + β2 * revenue + β3 * expenditure + β4 * university graduates + β5 * inflation + ε
unemployment = registered unemployment rate [%], share of registered unemployed in the civilian active population, i.e. excluding employees of budgetary units carrying out national defence and public security activities. (https://stat.gov.pl/obszary-tematyczne/rynek-pracy/bezrobocie-rejestrowane/stopa-bezrobocia-rejestrowanego-w-latach-1990-2025,4,1.html)
GDPpc = Gross Domestic Product per capita [zł], the end result of the production activities of resident production units (https://bdl.stat.gov.pl/bdl/dane/podgrup/tablica)
revenue = revenue of voivodships budgets per 1 inhabitant [zł] (https://bdl.stat.gov.pl/bdl/dane/podgrup/tablica)
expenditure = expenditure of voivodships budgets per 1 inhabitant [zł] (https://bdl.stat.gov.pl/bdl/dane/podgrup/tablica)
university_graduates = university graduates per 10,000 inhabitants (https://bdl.stat.gov.pl/bdl/dane/podgrup/tablica)
inflation = annual consumer price index, previous year = 100, (https://stat.gov.pl/obszary-tematyczne/ceny-handel/wskazniki-cen/wskazniki-cen-towarow-i-uslug-konsumpcyjnych-pot-inflacja-/roczne-wskazniki-cen-towarow-i-uslug-konsumpcyjnych/)
Analysis of the statistical significance of explanatory variables
To carry out the forecast, an initial linear model was developed as a basis for determining the subsequent modeling steps. A significance level of 0.05 was applied in both model construction and statistical hypothesis testing.
Linear regression model
With an assumed significance level of α =0.05, it should be concluded that the two explanatory variables university_graduates and inflation are not significant.
# model
model1 <- lm(unemployment ~ GDPpc + revenue + expenditure + university_graduates + inflation, data = df)
summary(model1) ##
## Call:
## lm(formula = unemployment ~ GDPpc + revenue + expenditure + university_graduates +
## inflation, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.89675 -0.39528 -0.02427 0.10495 1.15184
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.587e+01 1.023e+01 -2.530 0.03528 *
## GDPpc -4.341e-04 1.295e-04 -3.352 0.01005 *
## revenue 9.000e-02 2.566e-02 3.507 0.00799 **
## expenditure -6.163e-02 2.097e-02 -2.939 0.01874 *
## university_graduates 8.913e-02 4.423e-02 2.015 0.07866 .
## inflation 3.444e-01 1.604e-01 2.148 0.06402 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7197 on 8 degrees of freedom
## Multiple R-squared: 0.9715, Adjusted R-squared: 0.9536
## F-statistic: 54.48 on 5 and 8 DF, p-value: 5.776e-06Logarithmisation of variables
To improve model accuracy and reduce the p-values of the variables, all variables were log-transformed. Logarithmic models are often more precise and allow for the inclusion of more variables in the analysis.
variable_to_log <- c("unemployment", "GDPpc", "revenue", "expenditure", "university_graduates", "inflation")
df_log <- df
df_log[variable_to_log] <- lapply(df[variable_to_log], log)
model_log <- lm(unemployment ~ GDPpc + revenue + expenditure + university_graduates + inflation, data=df_log)
summary(model_log)##
## Call:
## lm(formula = unemployment ~ GDPpc + revenue + expenditure + university_graduates +
## inflation, data = df_log)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.07838 -0.02887 0.00470 0.03001 0.08701
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3151 3.2928 0.703 0.501951
## GDPpc -2.8509 0.5118 -5.570 0.000528 ***
## revenue 6.0347 1.0762 5.607 0.000506 ***
## expenditure -4.4188 0.8966 -4.928 0.001152 **
## university_graduates 1.0540 0.3489 3.021 0.016529 *
## inflation 3.4512 1.1342 3.043 0.015996 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05971 on 8 degrees of freedom
## Multiple R-squared: 0.9855, Adjusted R-squared: 0.9764
## F-statistic: 108.4 on 5 and 8 DF, p-value: 3.973e-07After log-transforming the variables, their p-values decreased. All explanatory variables proved to be statistically significant at the α = 0.05 significance level. Since both the dependent and independent variables are expressed in natural logarithms, the estimated coefficients can be interpreted as elasticities.
The coefficient for GDP per capita is -2.85 (p < 0.001), suggesting that a 1% increase in GDP per capita is associated with an approximate 2.85% decrease in the unemployment rate, ceteris paribus. This result is highly statistically significant and aligns with standard economic theory that links higher economic output with lower unemployment.
Similarly, public expenditure has a negative and significant effect on unemployment. The coefficient of −4.42 (p <0.05) implies that a 1% increase in voivodships spendings per 1 inhabitant leads to a 4.42% reduction in the unemployment rate, ceteris paribus.
Conversely, public revenue shows a positive elasticity of 6.03 (p < 0.001), indicating that a 1% increase in revenue of voivodships budgets per 1 inhabitant is associated with a 6.03% increase in unemployment rate, ceteris paribus.
The coefficient for university graduates is 1.05 (p <0,05), implying that a 1% rise in the number of university graduates per 10 000 inhabitants correlates with a 1.05% increase in unemployment rate, holding other factors constant. This unexpected finding might suggest a mismatch between higher education output and labor market demand.
Finally, inflation has a significant and positive coefficient of 3.45 (p = 0.0160), indicating that a 1% increase in the inflation rate is associated with a 3.45% increase in unemployment rate. This may reflect cost-push inflationary pressures within the economy.
The model demonstrates a high degree of explanatory power, with an R-squared of 0.9855 and an adjusted R-squared of 0.9764, indicating that approximately 98% of the variation in unemployment is explained by the included variables. The overall model is statistically significant (F-statistic = 108.4, p < 0.001).
Calculation of ex-ante errors
# utworzenie funkcji wyliczajacej blad RMSPE
fun_rmspe <- function(actual, predicted) {
rmse_value <- sqrt(mean((actual - predicted)^2))
rmspe_value <- (rmse_value / mean(actual)) * 100
}
# obliczanie bledow
MAE <- mean(abs(df_log$unemployment - model_log$fitted.values))
MSE <- mean((df_log$unemployment - model_log$fitted.values)^2)
RMSE <- sqrt(MSE)
MAPE <- mean(abs((df_log$unemployment - model_log$fitted.values) / df_log$unemployment)) * 100
rmspe_value <- fun_rmspe(df_log$unemployment, model_log$fitted.values)
# wyświetlenie bledow
cat("Mean Absolute Error (MAE):", MAE, "\n")## Mean Absolute Error (MAE): 0.03673092cat("Mean Squared Error (MSE):", MSE, "\n")## Mean Squared Error (MSE): 0.002037493cat("Root Mean Squared Error (RMSE):", RMSE, "\n")## Root Mean Squared Error (RMSE): 0.0451386cat("Mean Absolute Percentage Error (MAPE):", MAPE, "%\n")## Mean Absolute Percentage Error (MAPE): 1.83498 %cat("Root Mean Squared Percentage Error (RMSPE): ", rmspe_value, "\n")## Root Mean Squared Percentage Error (RMSPE): 2.163164The RMSPE for the predicted values of the dependent variable (unemployment) was approximately 2.16. This means that the predicted values differ from the actual ones by only 2.16% on average, indicating high model accuracy. The model predicts the dependent variable with a very low average percentage error, showing that the forecasted values are very close to the real ones.
Visualisation of logarithmic variables
## wizualizacja wykresów zmiennych
# utworzenie zbioru danych
df_ost <- subset(df_log, select = c(unemployment, GDPpc, revenue, expenditure, university_graduates, inflation))
lista_wykresow <- list()
# tworzenie wykresów dla zmiennej objaśnianej i zmiennych objaśniających
for (colname in colnames(df_ost)) {
nazwa_wykresu <- paste("wykres_", colname, sep = "")
# Tworzenie wykresu i zapisanie go do listy
assign(nazwa_wykresu,
ggplot(df_ost, aes_string(x = df_log$rok, y = colname)) +
geom_line() +
geom_point() +
labs(title = paste(colname, "in time"), x = "year", y = colname) +
theme_minimal()
)
lista_wykresow[[nazwa_wykresu]] <- get(nazwa_wykresu) #odanie utworzonego wykresu do listy z wykresami
}## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# wyświetlenie wszystkich wykresów z listy
for (wykres in lista_wykresow) {
print(wykres)
}
Statistical tests
Correlation
The correlation matrix shows that the unemployment variable is strongly negatively correlated with GDPpc (-0.89), public revenue (-0.67), public expenditure (-0.55) and inflation (-0.49), indicating that higher levels of GDPpc, public revenue and expenditure and inflation are associated with lower unemployment. The negative autocorrelation of these variables reflects classic macroeconomic relationships.
In addition, the high positive correlation of unemployment with the number of university graduates (0.95) is indicative of specific labour market phenomena, e.g. skills mismatches or delays in labour force absorption. This means that a higher number of university graduates is associated with a higher unemployment rate.
Strong positive correlations between economic variables (e.g. GDPpc with revenue and expenditure) indicate a high interdependence of these indicators and a potential risk of multicollinearity in the regression model.
library(corrplot)
# Macierz korelacji
cor_matrix <- cor(df_ost, use = "complete.obs")
# Pełna macierz kwadratowa
corrplot(cor_matrix,
method = "color",
type = "full",
col = colorRampPalette(c("#660fb8", "white", "#cf36c2"))(200),
addCoef.col = "white",
tl.col = "black",
tl.srt = 45,
number.cex = 0.7,
mar = c(0, 0, 2, 0)) 
Test for non-linearity - squares
The null hypothesis of the test for nonlinearity states that there is no nonlinearity present in the data that could be detected by adding squared terms to the model, while the alternative hypothesis assumes the presence of such nonlinearity. The obtained p-value of 0.21 exceeds the adopted significance level of 0.05, which means there is no basis to reject the null hypothesis. Therefore, the model does not show evidence of nonlinearity, and there is no justification for including squared terms in the model.
#test na nieliniowość (kwadraty)
model_log_test_kwadraty <- lm(unemployment ~ GDPpc + revenue + expenditure + university_graduates + inflation + I(GDPpc^2) + I(revenue^2) + I(expenditure^2) + I(university_graduates^2 + I(inflation)^2) , data = df_log)
anova(model_log, model_log_test_kwadraty)## Analysis of Variance Table
##
## Model 1: unemployment ~ GDPpc + revenue + expenditure + university_graduates +
## inflation
## Model 2: unemployment ~ GDPpc + revenue + expenditure + university_graduates +
## inflation + I(GDPpc^2) + I(revenue^2) + I(expenditure^2) +
## I(university_graduates^2 + I(inflation)^2)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 8 0.0285249
## 2 4 0.0083873 4 0.020138 2.401 0.2085Test for non-linearity (logarithms)
The null hypothesis of the test for nonlinearity assumes that there is no nonlinearity in the data that could be detected by adding logarithmic transformations of the variables, while the alternative hypothesis suggests the presence of such nonlinearity. The obtained p-value in the test is 0.002, which is significantly lower than the adopted significance level of 0.05. Therefore, the null hypothesis is rejected in favor of the alternative. This indicates the presence of nonlinearity, and it is justified to include logarithmic transformations of the variables in the model.
model_log_test_log <- lm(unemployment ~ log(GDPpc) + log(revenue) + log(expenditure) + log(university_graduates) + log(inflation), data = df_log)
anova(model_log, model_log_test_log)## Analysis of Variance Table
##
## Model 1: unemployment ~ GDPpc + revenue + expenditure + university_graduates +
## inflation
## Model 2: unemployment ~ log(GDPpc) + log(revenue) + log(expenditure) +
## log(university_graduates) + log(inflation)
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 8 0.028525
## 2 8 0.026455 0 0.0020703Durbin Watson test
The Durbin Watson test is used to detect the first-order autocorrelation of residuals. Its null hypothesis states that autocorrelation of the residuals does not occur, while the alternative hypothesis indicates its presence. The p-value is 0.48, which is greater than the significance level of 0.05. Therefore, there is no basis to reject the null hypothesis, suggesting that the residuals are not autocorrelated
dwtest(model_log)##
## Durbin-Watson test
##
## data: model_log
## DW = 2.5679, p-value = 0.4758
## alternative hypothesis: true autocorrelation is greater than 0Breusch-Godfrey test
The Breusch-Godfrey test is used to detect the serial autocorrelation of the residuals. The null hypothesis states that there is no autocorrelation, while the alternative hypothesis indicates its presence. The p-value is 0.09, which is greater than the significance level of 0.05. Therefore, there is no sufficient evidence to reject the null hypothesis, suggesting that the residuals are not autocorrelated.
bgtest(model_log)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_log
## LM test = 2.8824, df = 1, p-value = 0.08955Breusch-Pagan test
The Breusch-Pagan test is used to detect whether heteroskedasticity occurs in the model. It is an undesirable situation in which the variance of the residuals of the model are not constant. The null hypothesis assumes homoskedasticity, meaning the variance of the error term is constant, while the alternative hypothesis indicates heteroskedasticity. The p-value of 0.08 exceeds the significance level of 0.05, providing no sufficient evidence to reject the null hypothesis. Therefore, it can be concluded that the model does not exhibit heteroskedasticity, and the assumption of constant variance of the error term holds.
bptest(model_log)##
## studentized Breusch-Pagan test
##
## data: model_log
## BP = 9.7198, df = 5, p-value = 0.08358The RESET (Regression Equation Specification Error Test) test
The RESET test is used to detect specification errors in a linear regression model. It assesses whether the assumed functional form of the model is correct. The null hypothesis assumes that the model is correctly specified, while the alternative suggests misspecification. The p-value of 0.25 exceeds the significance level of 0.05, providing no basis to reject the null hypothesis. Thus, the model appears to be correctly specified.
resettest(model_log)##
## RESET test
##
## data: model_log
## RESET = 1.7843, df1 = 2, df2 = 6, p-value = 0.2466Shapiro-Wilk test
The Shapiro-Wilk test is used to assess whether the residuals of the model follow a normal distribution. The null hypothesis assumes that the residuals are normally distributed, while the alternative hypothesis suggests a deviation from normality. The p-value is 0.97, which is significantly higher than the significance level of 0.05. Therefore, there is no basis to reject the null hypothesis, indicating that the residuals are normally distributed.
shapiro.test(resid(model_log))##
## Shapiro-Wilk normality test
##
## data: resid(model_log)
## W = 0.97868, p-value = 0.9664Forecasting of Explanatory Variables
Indeks = 2:15
df_log <- df_log %>%
mutate(Indeks) %>%
select(Indeks, everything())# tworzenie zmiennych
log_GDPpc <- df_log$GDPpc
log_revenue <- df_log$revenue
log_expenditure <- df_log$expenditure
log_university_graduates <- df_log$university_graduates
log_inflation <- df_log$inflation
# zamienienie tych zmiennych na szeregi czasowe
log_GDPpc.ts <- ts(log_GDPpc, frequency=1, start=c(2010))
log_revenue.ts <- ts(log_revenue, frequency=1, start=c(2010))
log_expenditure.ts <- ts(log_expenditure, frequency=1, start=c(2010))
log_university_graduates.ts <- ts(log_university_graduates, frequency=1, start=c(2010))
log_inflation.ts <- ts(log_inflation, frequency=1, start=c(2010))
# utworzenie zbioru danych z utworzonych szeregow czasowych
df_ts <- cbind(log_GDPpc.ts, log_revenue.ts, log_expenditure.ts, log_university_graduates.ts, log_inflation.ts)
okres_prognozy <- data.frame(Indeks = seq(20,20))Forecasting- Linear Trend
## PROGNOZOWANIE - TREND LINIOWY
# Upewnij się, że indeks jest od 1 do 14 dla 2010–2023
df_log$Indeks <- 1:nrow(df_log)
# Nazwy kolumn w df_ts
colnames(df_ts) <- c("log_GDPpc", "log_revenue", "log_expenditure", "log_university_graduates", "log_inflation")
# Dane prognozowane – indeksy 15:21 (czyli 2024–2030)
okres_prognozy <- data.frame(Indeks = 15:21)
# Funkcja RMSPE
fun_rmspe <- function(actual, predicted) {
sqrt(mean((actual - predicted)^2)) / mean(actual) * 100
}
# Listy na wyniki
rmspe_values_lm <- list()
lista_wykresow_lm <- list()
for (var in colnames(df_ts)) {
cat("Model:", var, "\n")
df_model <- data.frame(
Indeks = df_log$Indeks,
Y = df_ts[, var]
)
# Model
model_progn <- lm(Y ~ Indeks, data = df_model)
summary(model_progn)
# Prognoza – zwracana jako matrix
prognoza <- predict(model_progn, newdata = okres_prognozy, interval = "confidence")
# RMSPE
fitted_vals <- model_progn$fitted.values
actual_vals <- df_model$Y
rmspe <- fun_rmspe(actual_vals, fitted_vals)
rmspe_values_lm[[var]] <- round(rmspe, 2)
# Funkcja MASE
fun_mase <- function(actual, predicted, training) {
n <- length(training)
scale <- mean(abs(diff(training)))
mean(abs(actual - predicted)) / scale
}
# Tworzenie wektora lat i danych do wykresu
Rok <- c(2010:2023, 2024:2030)
Actual <- c(actual_vals, rep(NA, 7))
Predicted <- c(fitted_vals, prognoza[, "fit"])
Lower <- c(rep(NA, length(fitted_vals)), prognoza[, "lwr"])
Upper <- c(rep(NA, length(fitted_vals)), prognoza[, "upr"])
wykres_df <- data.frame(Rok, Actual, Predicted, Lower, Upper)
# Tworzenie wykresu
wykres <- ggplot(wykres_df, aes(x = Rok)) +
geom_line(aes(y = Actual), color = "pink", size = 1) +
geom_line(aes(y = Predicted), color = "purple", linetype = "dashed", size = 1) +
geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "purple", alpha = 0.2) +
labs(title = paste("Linear trend forecast for:", var),
y = var, x = "year") +
theme_minimal()
lista_wykresow_lm[[var]] <- wykres
cat("RMSPE:", round(rmspe, 2), "\n--------------------------\n")
}## Model: log_GDPpc## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## RMSPE: 0.54
## --------------------------
## Model: log_revenue
## RMSPE: 2.2
## --------------------------
## Model: log_expenditure
## RMSPE: 2.59
## --------------------------
## Model: log_university_graduates
## RMSPE: 1.07
## --------------------------
## Model: log_inflation
## RMSPE: 0.69
## --------------------------results <- list()
for (var in colnames(df_ts)) {
cat("Model dfor the variable:", var, "\n")
# Dane do modelu
df_model <- data.frame(Indeks = df_log$Indeks, Y = df_ts[, var])
# Model liniowy
model_progn <- lm(Y ~ Indeks, data = df_model)
print(summary(model_progn))
# Testy statystyczne
cat("Test Durbin-Watson:\n")
print(dwtest(model_progn))
cat("Test Breusch-Godfrey:\n")
print(bgtest(model_progn))
cat("Test Breusch-Pagan:\n")
print(bptest(model_progn))
cat("Test RESET:\n")
print(resettest(model_progn))
# Prognoza 2024-2030
prognoza <- predict(model_progn, newdata = okres_prognozy, interval = "confidence")
# Miary dokładności
fitted_vals <- model_progn$fitted.values
actual_vals <- df_model$Y
ME <- mean(actual_vals - fitted_vals)
RMSE <- sqrt(mean((actual_vals - fitted_vals)^2))
MAE <- mean(abs(actual_vals - fitted_vals))
MAPE <- mean(abs((actual_vals - fitted_vals) / actual_vals)) * 100
MASE <- fun_mase(actual_vals, fitted_vals, actual_vals)
RMSPE <- fun_rmspe(actual_vals, fitted_vals)
cat("Accuracy metrics on the training set:\n")
cat("ME:", round(ME, 4), "\n")
cat("RMSE:", round(RMSE, 4), "\n")
cat("MAE:", round(MAE, 4), "\n")
cat("MAPE (%):", round(MAPE, 2), "\n")
cat("MASE:", round(MASE, 4), "\n")
cat("RMSPE (%):", round(RMSPE, 2), "\n")
cat("-----------------------------\n\n")
# Prognoza
prognoza_df <- data.frame(
Rok = 2024:2030,
Prognoza = prognoza[, "fit"],
Dolny_przedzial = prognoza[, "lwr"],
Gorny_przedzial = prognoza[, "upr"]
)
cat("Forecast for the variable", var, "for the years 2024-2030:\n")
print(prognoza_df)
cat("\n==============================\n\n")
# Zapis do listy wyników
results[[var]] <- list(
model_summary = summary(model_progn),
dwtest = dwtest(model_progn),
bgtest = bgtest(model_progn),
bptest = bptest(model_progn),
resettest = resettest(model_progn),
forecast = prognoza_df,
accuracy_measures = data.frame(ME = ME, RMSE = RMSE, MAE = MAE, MAPE = MAPE, MASE = MASE, RMSPE = RMSPE)
)
}## Model dfor the variable: log_GDPpc
##
## Call:
## lm(formula = Y ~ Indeks, data = df_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06990 -0.04996 -0.01860 0.05304 0.11002
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.405214 0.035909 289.76 < 2e-16 ***
## Indeks 0.064052 0.004217 15.19 3.37e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06361 on 12 degrees of freedom
## Multiple R-squared: 0.9506, Adjusted R-squared: 0.9464
## F-statistic: 230.7 on 1 and 12 DF, p-value: 3.371e-09
##
## Test Durbin-Watson:
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.30229, p-value = 6.226e-08
## alternative hypothesis: true autocorrelation is greater than 0
##
## Test Breusch-Godfrey:
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 9.9889, df = 1, p-value = 0.001575
##
## Test Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 2.4562, df = 1, p-value = 0.1171
##
## Test RESET:
##
## RESET test
##
## data: model_progn
## RESET = 92.666, df1 = 2, df2 = 10, p-value = 3.517e-07
##
## Accuracy metrics on the training set:
## ME: 0
## RMSE: 0.0589
## MAE: 0.0528
## MAPE (%): 0.48
## MASE: 0.7744
## RMSPE (%): 0.54
## -----------------------------
##
## Forecast for the variable log_GDPpc for the years 2024-2030:
## Rok Prognoza Dolny_przedzial Gorny_przedzial
## 1 2024 11.36599 11.28775 11.44423
## 2 2025 11.43004 11.34360 11.51649
## 3 2026 11.49410 11.39927 11.58892
## 4 2027 11.55815 11.45480 11.66150
## 5 2028 11.62220 11.51023 11.73417
## 6 2029 11.68625 11.56557 11.80694
## 7 2030 11.75030 11.62084 11.87976
##
## ==============================
##
## Model dfor the variable: log_revenue
##
## Call:
## lm(formula = Y ~ Indeks, data = df_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26660 -0.06963 0.03558 0.06889 0.27163
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.780754 0.082544 70.032 < 2e-16 ***
## Indeks 0.049707 0.009694 5.127 0.00025 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1462 on 12 degrees of freedom
## Multiple R-squared: 0.6866, Adjusted R-squared: 0.6605
## F-statistic: 26.29 on 1 and 12 DF, p-value: 0.0002502
##
## Test Durbin-Watson:
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.62421, p-value = 0.0001769
## alternative hypothesis: true autocorrelation is greater than 0
##
## Test Breusch-Godfrey:
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 6.4265, df = 1, p-value = 0.01124
##
## Test Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 1.1227, df = 1, p-value = 0.2893
##
## Test RESET:
##
## RESET test
##
## data: model_progn
## RESET = 13.493, df1 = 2, df2 = 10, p-value = 0.001445
##
## Accuracy metrics on the training set:
## ME: 0
## RMSE: 0.1354
## MAE: 0.1053
## MAPE (%): 1.72
## MASE: 0.9864
## RMSPE (%): 2.2
## -----------------------------
##
## Forecast for the variable log_revenue for the years 2024-2030:
## Rok Prognoza Dolny_przedzial Gorny_przedzial
## 1 2024 6.526351 6.346503 6.706199
## 2 2025 6.576058 6.377353 6.774763
## 3 2026 6.625764 6.407787 6.843742
## 4 2027 6.675471 6.437906 6.913036
## 5 2028 6.725177 6.467782 6.982572
## 6 2029 6.774884 6.497468 7.052300
## 7 2030 6.824590 6.527001 7.122180
##
## ==============================
##
## Model dfor the variable: log_expenditure
##
## Call:
## lm(formula = Y ~ Indeks, data = df_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.33283 -0.06644 0.03123 0.08919 0.32975
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.84803 0.09710 60.225 2.9e-16 ***
## Indeks 0.03984 0.01140 3.493 0.00444 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.172 on 12 degrees of freedom
## Multiple R-squared: 0.5042, Adjusted R-squared: 0.4629
## F-statistic: 12.2 on 1 and 12 DF, p-value: 0.004435
##
## Test Durbin-Watson:
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.67453, p-value = 0.00034
## alternative hypothesis: true autocorrelation is greater than 0
##
## Test Breusch-Godfrey:
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 5.8105, df = 1, p-value = 0.01593
##
## Test Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 1.0304, df = 1, p-value = 0.3101
##
## Test RESET:
##
## RESET test
##
## data: model_progn
## RESET = 13.719, df1 = 2, df2 = 10, p-value = 0.00136
##
## Accuracy metrics on the training set:
## ME: 0
## RMSE: 0.1593
## MAE: 0.1219
## MAPE (%): 1.99
## MASE: 1.0155
## RMSPE (%): 2.59
## -----------------------------
##
## Forecast for the variable log_expenditure for the years 2024-2030:
## Rok Prognoza Dolny_przedzial Gorny_przedzial
## 1 2024 6.445634 6.234063 6.657204
## 2 2025 6.485474 6.251720 6.719227
## 3 2026 6.525314 6.268889 6.781739
## 4 2027 6.565154 6.285687 6.844622
## 5 2028 6.604994 6.302199 6.907790
## 6 2029 6.644835 6.318487 6.971182
## 7 2030 6.684675 6.334595 7.034754
##
## ==============================
##
## Model dfor the variable: log_university_graduates
##
## Call:
## lm(formula = Y ~ Indeks, data = df_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.069467 -0.049763 0.002907 0.034895 0.080578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.916124 0.029832 164.79 < 2e-16 ***
## Indeks -0.045714 0.003504 -13.05 1.89e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05284 on 12 degrees of freedom
## Multiple R-squared: 0.9342, Adjusted R-squared: 0.9287
## F-statistic: 170.2 on 1 and 12 DF, p-value: 1.893e-08
##
## Test Durbin-Watson:
##
## Durbin-Watson test
##
## data: model_progn
## DW = 1.3935, p-value = 0.06061
## alternative hypothesis: true autocorrelation is greater than 0
##
## Test Breusch-Godfrey:
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 0.45643, df = 1, p-value = 0.4993
##
## Test Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 1.023, df = 1, p-value = 0.3118
##
## Test RESET:
##
## RESET test
##
## data: model_progn
## RESET = 7.8781, df1 = 2, df2 = 10, p-value = 0.008822
##
## Accuracy metrics on the training set:
## ME: 0
## RMSE: 0.0489
## MAE: 0.0432
## MAPE (%): 0.95
## MASE: 0.7913
## RMSPE (%): 1.07
## -----------------------------
##
## Forecast for the variable log_university_graduates for the years 2024-2030:
## Rok Prognoza Dolny_przedzial Gorny_przedzial
## 1 2024 4.230417 4.165419 4.295415
## 2 2025 4.184704 4.112891 4.256516
## 3 2026 4.138990 4.060212 4.217768
## 4 2027 4.093276 4.007419 4.179133
## 5 2028 4.047562 3.954539 4.140586
## 6 2029 4.001848 3.901589 4.102108
## 7 2030 3.956135 3.848585 4.063685
##
## ==============================
##
## Model dfor the variable: log_inflation
##
## Call:
## lm(formula = Y ~ Indeks, data = df_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.03760 -0.02469 -0.01391 0.02809 0.06875
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.596867 0.019602 234.509 <2e-16 ***
## Indeks 0.005698 0.002302 2.475 0.0292 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03472 on 12 degrees of freedom
## Multiple R-squared: 0.338, Adjusted R-squared: 0.2828
## F-statistic: 6.127 on 1 and 12 DF, p-value: 0.02921
##
## Test Durbin-Watson:
##
## Durbin-Watson test
##
## data: model_progn
## DW = 0.67131, p-value = 0.0003267
## alternative hypothesis: true autocorrelation is greater than 0
##
## Test Breusch-Godfrey:
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: model_progn
## LM test = 5.5599, df = 1, p-value = 0.01838
##
## Test Breusch-Pagan:
##
## studentized Breusch-Pagan test
##
## data: model_progn
## BP = 1.2547, df = 1, p-value = 0.2627
##
## Test RESET:
##
## RESET test
##
## data: model_progn
## RESET = 15.548, df1 = 2, df2 = 10, p-value = 0.0008531
##
## Accuracy metrics on the training set:
## ME: 0
## RMSE: 0.0321
## MAE: 0.0286
## MAPE (%): 0.61
## MASE: 1.5129
## RMSPE (%): 0.69
## -----------------------------
##
## Forecast for the variable log_inflation for the years 2024-2030:
## Rok Prognoza Dolny_przedzial Gorny_przedzial
## 1 2024 4.682343 4.639633 4.725052
## 2 2025 4.688041 4.640854 4.735228
## 3 2026 4.693739 4.641975 4.745503
## 4 2027 4.699438 4.643022 4.755853
## 5 2028 4.705136 4.644011 4.766261
## 6 2029 4.710834 4.644955 4.776714
## 7 2030 4.716533 4.645863 4.787203
##
## ==============================# Wykresy
for (wykres in lista_wykresow_lm) {
print(wykres)
}## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
Forecasting - ARIMA
# Inicjalizacja list
rmspe_values_arima <- list()
lista_wykresow_arima <- list()
# Inicjalizacja ramki do zbierania prognoz
prognozy_arima_df <- data.frame(rok = 2010:2030)
# Pętla po zmiennych objaśniających
for (var in colnames(df_log[, -c(1, 2)])) {
cat("\nModel\n")
cat(paste0("[1] \"", var, "\"\n\n"))
# dopasowanie modelu ARIMA
arima_model <- auto.arima(df_log[[var]])
print(summary(arima_model))
# prognoza
arima_prognoza <- forecast(arima_model, h = 7)
print(arima_prognoza)
# RMSPE
rmspe_value <- fun_rmspe(df_log[[var]], arima_prognoza$fitted)
cat("RMSPE: ", round(rmspe_value, 6), "\n")
cat("-------------------------------------------------------------------------------\n")
rmspe_values_arima[[var]] <- paste("RMSPE (ARIMA):", round(rmspe_value, 6))
# Testy reszt
cat("[1] \"Tests for residuals\"\n")
cat("ADF test:\n"); print(adf.test(residuals(arima_model)))
cat("Box test:\n"); print(Box.test(residuals(arima_model), lag = 10, type = "Ljung-Box"))
# Dane do wykresu
wykres_dane_arima <- data.frame(
rok = c(df_log$rok, 2024:2030),
actual = c(df_log[[var]], rep(NA, 7)),
predicted = c(arima_prognoza$fitted, arima_prognoza$mean)
)
wykres_prognozowane_dane_arima <- data.frame(
Time = 2024:2030,
Forecast = as.numeric(arima_prognoza$mean),
Lo80 = arima_prognoza$lower[, "80%"],
Hi80 = arima_prognoza$upper[, "80%"],
Lo95 = arima_prognoza$lower[, "95%"],
Hi95 = arima_prognoza$upper[, "95%"]
)
# Wykres
wykres <- ggplot(wykres_dane_arima, aes(x = rok)) +
geom_line(aes(y = actual), color = "pink", size = 1) +
geom_line(aes(y = predicted), color = "purple", linetype = "dashed", size = 1) +
geom_ribbon(data = wykres_prognozowane_dane_arima,
aes(x = Time, ymin = Lo95, ymax = Hi95), fill = "purple", alpha = 0.2) +
geom_ribbon(data = wykres_prognozowane_dane_arima,
aes(x = Time, ymin = Lo80, ymax = Hi80), fill = "purple", alpha = 0.4) +
labs(title = paste("ARIMA Forecast:", var),
x = "year", y = var) +
theme_minimal()
print(wykres)
# Zapis wykresu do listy
nazwa_wykresu <- paste("wykres_arima_", var, sep = "")
assign(nazwa_wykresu, wykres)
lista_wykresow_arima[[nazwa_wykresu]] <- wykres
# Zbieranie fitted + forecast do ramki prognozy
prognoza_pelna <- c(as.numeric(arima_prognoza$fitted), as.numeric(arima_prognoza$mean))
stopifnot(length(prognoza_pelna) == length(2010:2030))
prognozy_arima_df[[var]] <- prognoza_pelna
}##
## Model
## [1] "unemployment"
##
## Series: df_log[[var]]
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## -0.0683
## s.e. 0.0304
##
## sigma^2 = 0.01303: log likelihood = 10.29
## AIC=-16.58 AICc=-15.38 BIC=-15.45
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.000184717 0.1056651 0.08471971 -0.1122364 4.142828 0.7716067
## ACF1
## Training set 0.2708364
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 15 1.560898 1.4146325 1.707163 1.3372044 1.784591
## 16 1.492555 1.2857047 1.699405 1.1762048 1.808905
## 17 1.424212 1.1708734 1.677551 1.0367639 1.811661
## 18 1.355869 1.0633390 1.648400 0.9084827 1.803256
## 19 1.287527 0.9604677 1.614586 0.7873331 1.787720
## 20 1.219184 0.8609087 1.577459 0.6712493 1.767119
## 21 1.150841 0.7638597 1.537823 0.5590041 1.742678
## RMSPE: 5.06376
## -------------------------------------------------------------------------------
## [1] "Tests for residuals"
## ADF test:
##
## Augmented Dickey-Fuller Test
##
## data: residuals(arima_model)
## Dickey-Fuller = -1.8888, Lag order = 2, p-value = 0.6147
## alternative hypothesis: stationary
##
## Box test:
##
## Box-Ljung test
##
## data: residuals(arima_model)
## X-squared = 10.577, df = 10, p-value = 0.3914## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
##
## Model
## [1] "GDPpc"
##
## Series: df_log[[var]]
## ARIMA(0,2,0)
##
## sigma^2 = 0.0009043: log likelihood = 25.14
## AIC=-48.29 AICc=-47.89 BIC=-47.8
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.000546014 0.02784098 0.02233826 0.004162227 0.204867 0.3274605
## ACF1
## Training set -0.04080073
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 15 11.50847 11.46993 11.54701 11.44953 11.56741
## 16 11.60499 11.51882 11.69117 11.47320 11.73678
## 17 11.70151 11.55731 11.84571 11.48098 11.92204
## 18 11.79803 11.58694 12.00911 11.47520 12.12085
## 19 11.89454 11.60874 12.18035 11.45744 12.33165
## 20 11.99106 11.62343 12.35870 11.42882 12.55331
## 21 12.08758 11.63159 12.54357 11.39020 12.78496
## RMSPE: 0.25576
## -------------------------------------------------------------------------------
## [1] "Tests for residuals"
## ADF test:
##
## Augmented Dickey-Fuller Test
##
## data: residuals(arima_model)
## Dickey-Fuller = -2.5188, Lag order = 2, p-value = 0.3747
## alternative hypothesis: stationary
##
## Box test:
##
## Box-Ljung test
##
## data: residuals(arima_model)
## X-squared = 14.088, df = 10, p-value = 0.169## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
##
## Model
## [1] "revenue"
##
## Series: df_log[[var]]
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## 0.0650
## s.e. 0.0305
##
## sigma^2 = 0.01309: log likelihood = 10.26
## AIC=-16.51 AICc=-15.31 BIC=-15.38
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0004170079 0.1059447 0.06661939 -0.03600442 1.080687 0.6238316
## ACF1
## Training set 0.198698
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 15 6.813291 6.666638 6.959943 6.589005 7.037576
## 16 6.878302 6.670905 7.085700 6.561115 7.195490
## 17 6.943314 6.689305 7.197323 6.554841 7.331788
## 18 7.008326 6.715022 7.301630 6.559756 7.456896
## 19 7.073338 6.745414 7.401262 6.571821 7.574855
## 20 7.138350 6.779127 7.497573 6.588966 7.687734
## 21 7.203362 6.815356 7.591367 6.609959 7.796764
## RMSPE: 1.721683
## -------------------------------------------------------------------------------
## [1] "Tests for residuals"
## ADF test:
##
## Augmented Dickey-Fuller Test
##
## data: residuals(arima_model)
## Dickey-Fuller = -1.5156, Lag order = 2, p-value = 0.7569
## alternative hypothesis: stationary
##
## Box test:
##
## Box-Ljung test
##
## data: residuals(arima_model)
## X-squared = 7.4372, df = 10, p-value = 0.6836## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
##
## Model
## [1] "expenditure"
##
## Series: df_log[[var]]
## ARIMA(0,1,0)
##
## sigma^2 = 0.02146: log likelihood = 6.52
## AIC=-11.05 AICc=-10.68 BIC=-10.48
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05432959 0.1411756 0.1119279 0.829781 1.813399 0.9321292
## ACF1
## Training set 0.1699349
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 15 6.735542 6.547789 6.923296 6.448398 7.022687
## 16 6.735542 6.470019 7.001066 6.329459 7.141626
## 17 6.735542 6.410344 7.060741 6.238194 7.232891
## 18 6.735542 6.360035 7.111050 6.161254 7.309831
## 19 6.735542 6.315713 7.155372 6.093468 7.377617
## 20 6.735542 6.275642 7.195443 6.032186 7.438899
## 21 6.735542 6.238793 7.232292 5.975830 7.495255
## RMSPE: 2.296722
## -------------------------------------------------------------------------------
## [1] "Tests for residuals"
## ADF test:
##
## Augmented Dickey-Fuller Test
##
## data: residuals(arima_model)
## Dickey-Fuller = -1.8063, Lag order = 2, p-value = 0.6462
## alternative hypothesis: stationary
##
## Box test:
##
## Box-Ljung test
##
## data: residuals(arima_model)
## X-squared = 8.6662, df = 10, p-value = 0.564## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
##
## Model
## [1] "university_graduates"
##
## Series: df_log[[var]]
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## -0.0350
## s.e. 0.0164
##
## sigma^2 = 0.00377: log likelihood = 18.35
## AIC=-32.7 AICc=-31.5 BIC=-31.57
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0003462299 0.05684441 0.04441252 0.001872576 0.9751763 0.8140261
## ACF1
## Training set -0.2633499
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 15 4.321672 4.242986 4.400358 4.201332 4.442012
## 16 4.286636 4.175357 4.397914 4.116449 4.456822
## 17 4.251599 4.115311 4.387887 4.043164 4.460034
## 18 4.216563 4.059191 4.373934 3.975883 4.457242
## 19 4.181526 4.005579 4.357473 3.912438 4.450614
## 20 4.146489 3.953749 4.339230 3.851718 4.441260
## 21 4.111453 3.903269 4.319636 3.793064 4.429842
## RMSPE: 1.242971
## -------------------------------------------------------------------------------
## [1] "Tests for residuals"
## ADF test:
##
## Augmented Dickey-Fuller Test
##
## data: residuals(arima_model)
## Dickey-Fuller = -1.971, Lag order = 2, p-value = 0.5834
## alternative hypothesis: stationary
##
## Box test:
##
## Box-Ljung test
##
## data: residuals(arima_model)
## X-squared = 7.8638, df = 10, p-value = 0.6421## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
##
## Model
## [1] "inflation"
##
## Series: df_log[[var]]
## ARIMA(0,1,0)
##
## sigma^2 = 0.0007885: log likelihood = 28.01
## AIC=-54.03 AICc=-53.66 BIC=-53.46
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.006208588 0.02705825 0.01787215 0.1312967 0.3823878 0.9460813
## ACF1
## Training set -0.09776298
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 15 4.713127 4.677142 4.749113 4.658092 4.768162
## 16 4.713127 4.662236 4.764019 4.635296 4.790959
## 17 4.713127 4.650799 4.775456 4.617804 4.808451
## 18 4.713127 4.641156 4.785098 4.603057 4.823198
## 19 4.713127 4.632661 4.793593 4.590065 4.836190
## 20 4.713127 4.624981 4.801274 4.578319 4.847935
## 21 4.713127 4.617919 4.808336 4.567518 4.858737
## RMSPE: 0.583202
## -------------------------------------------------------------------------------
## [1] "Tests for residuals"
## ADF test:
##
## Augmented Dickey-Fuller Test
##
## data: residuals(arima_model)
## Dickey-Fuller = -2.94, Lag order = 2, p-value = 0.2143
## alternative hypothesis: stationary
##
## Box test:
##
## Box-Ljung test
##
## data: residuals(arima_model)
## X-squared = 5.2723, df = 10, p-value = 0.8723## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
Forecasting - UNEMPLOYMENT
library(readxl)
library(tibble)
library(dplyr)
library(knitr)
# Wczytaj dane z Excela
df_excel <- read_excel("C:/Users/zuzal/OneDrive/Desktop/predicted vals.xlsx", sheet = "Arkusz2")
prognozy_final <- df_excel %>%
rename(Rok = rok, Zmienna = zmienna, Trend = trend, ARIMA = ARIMA) %>%
mutate(
Trend = round(Trend, 2),
ARIMA = round(ARIMA, 2)
)
kable(prognozy_final, caption = "Comparison of forecasts by method for all variebles")| Rok | Zmienna | Trend | ARIMA |
|---|---|---|---|
| 2024 | log_GDPpc | 11.37 | 11.51 |
| 2025 | log_GDPpc | 11.43 | 11.60 |
| 2026 | log_GDPpc | 11.49 | 11.70 |
| 2027 | log_GDPpc | 11.56 | 11.80 |
| 2028 | log_GDPpc | 11.62 | 11.89 |
| 2029 | log_GDPpc | 11.69 | 11.99 |
| 2030 | log_GDPpc | 11.75 | 12.09 |
| 2024 | log_revenue | 6.53 | 6.81 |
| 2025 | log_revenue | 6.58 | 6.88 |
| 2026 | log_revenue | 6.63 | 6.94 |
| 2027 | log_revenue | 6.68 | 7.01 |
| 2028 | log_revenue | 6.73 | 7.07 |
| 2029 | log_revenue | 6.77 | 7.14 |
| 2030 | log_revenue | 6.82 | 7.20 |
| 2024 | log_expenditure | 6.45 | 6.74 |
| 2025 | log_expenditure | 6.49 | 6.74 |
| 2026 | log_expenditure | 6.53 | 6.74 |
| 2027 | log_expenditure | 6.57 | 6.74 |
| 2028 | log_expenditure | 6.60 | 6.74 |
| 2029 | log_expenditure | 6.64 | 6.74 |
| 2030 | log_expenditure | 6.68 | 6.74 |
| 2024 | log_university_graduates | 4.23 | 4.32 |
| 2025 | log_university_graduates | 4.18 | 4.29 |
| 2026 | log_university_graduates | 4.14 | 4.25 |
| 2027 | log_university_graduates | 4.09 | 4.22 |
| 2028 | log_university_graduates | 4.05 | 4.18 |
| 2029 | log_university_graduates | 4.00 | 4.15 |
| 2030 | log_university_graduates | 3.96 | 4.11 |
| 2024 | log_inflation | 4.68 | 4.71 |
| 2025 | log_inflation | 4.69 | 4.71 |
| 2026 | log_inflation | 4.69 | 4.71 |
| 2027 | log_inflation | 4.70 | 4.71 |
| 2028 | log_inflation | 4.71 | 4.71 |
| 2029 | log_inflation | 4.71 | 4.71 |
| 2030 | log_inflation | 4.72 | 4.71 |
Forecasting unemployment- Linear trend
### PROGNOZOWANIE unemployment
# Dane do analizy
ts_unemployment <- ts(df_log$unemployment, start = 2010, frequency = 1)
future_years <- 2024:2030
rok_all <- c(2010:2023, future_years)
Indeks <- 1:length(rok_all)
h <- length(future_years)
# MODEL TREND LINIOWY
df_trend <- data.frame(Indeks = 1:length(ts_unemployment), Y = ts_unemployment)
model_trend <- lm(Y ~ Indeks, data = df_trend)
# prognoza
forecast_trend <- predict(model_trend, newdata = data.frame(Indeks = 15:21), interval = "confidence")
# dane do wykresu
df_trend_plot <- data.frame(
Rok = rok_all,
Actual = c(ts_unemployment, rep(NA, h)),
Predicted = c(model_trend$fitted.values, forecast_trend[, "fit"]),
Lower = c(rep(NA, length(ts_unemployment)), forecast_trend[, "lwr"]),
Upper = c(rep(NA, length(ts_unemployment)), forecast_trend[, "upr"])
)
ggplot(df_trend_plot, aes(x = Rok)) +
geom_line(aes(y = Actual), color = "pink") +
geom_line(aes(y = Predicted), color = "purple", linetype = "dashed") +
geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "purple", alpha = 0.2) +
labs(title = "Linear Trend Forecast – unemployment (log)", y = "log-unemployment", x = "year") +
theme_minimal()## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
# DIAGNOSTYKA MODEL_TREND
# Reszty
res_trend <- residuals(model_trend)
# Durbin-Watson Test
dwtest(model_trend)##
## Durbin-Watson test
##
## data: model_trend
## DW = 0.70467, p-value = 0.0004869
## alternative hypothesis: true autocorrelation is greater than 0# Breusch-Godfrey Test (2 rzędy opóźnień)
bgtest(model_trend, order = 2)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: model_trend
## LM test = 6.8732, df = 2, p-value = 0.03217# Breusch-Pagan Test
bptest(model_trend)##
## studentized Breusch-Pagan test
##
## data: model_trend
## BP = 0.089291, df = 1, p-value = 0.7651# RESET Test
resettest(model_trend)##
## RESET test
##
## data: model_trend
## RESET = 6.0693, df1 = 2, df2 = 10, p-value = 0.0188# MIARY DOKŁADNOŚCI – MODEL_TREND
# Prognoza vs rzeczywistość
actual_trend <- ts_unemployment
fitted_trend <- model_trend$fitted.values
me_trend <- mean(fitted_trend - actual_trend)
rmse_trend <- sqrt(mean((fitted_trend - actual_trend)^2))
mae_trend <- mean(abs(fitted_trend - actual_trend))
mape_trend <- mean(abs((fitted_trend - actual_trend) / actual_trend)) * 100
mase_trend <- mean(abs(fitted_trend - actual_trend)) / mean(abs(diff(actual_trend)))
rmspe_trend <- sqrt(mean(((fitted_trend - actual_trend) / actual_trend)^2)) * 100
cat("\n== Trend Model Accuracy ==\n")##
## == Trend Model Accuracy ==cat("ME:", me_trend, "\nRMSE:", rmse_trend, "\nMAE:", mae_trend, "\nMAPE:", mape_trend,
"\nMASE:", mase_trend, "\nRMSPE:", rmspe_trend, "\n")## ME: 1.586372e-17
## RMSE: 0.1277312
## MAE: 0.1110095
## MAPE: 5.517675
## MASE: 1.011048
## RMSPE: 6.529023# Prognozy – log(unemployment)
prognozy_log <- forecast_trend[, "fit"]
# Prognozy – odlogarytmowane
prognozy_real <- exp(prognozy_log)Forecasting unemployment - ARIMA
library(FinTS)## Warning: pakiet 'FinTS' został zbudowany w wersji R 4.4.3##
## Dołączanie pakietu: 'FinTS'## Następujący obiekt został zakryty z 'package:forecast':
##
## Acf# ARIMA
model_arima <- auto.arima(ts_unemployment)
forecast_arima <- forecast(model_arima, h = h)
# wykres
df_arima <- data.frame(
Rok = rok_all,
Actual = c(ts_unemployment, rep(NA, h)),
Predicted = c(rep(NA, length(ts_unemployment)), forecast_arima$mean),
Lower = c(rep(NA, length(ts_unemployment)), forecast_arima$lower[,2]),
Upper = c(rep(NA, length(ts_unemployment)), forecast_arima$upper[,2])
)
df_arima <- data.frame(
Rok = rok_all,
Actual = c(ts_unemployment, rep(NA, h)),
Predicted = c(fitted(model_arima), forecast_arima$mean),
Lower = c(rep(NA, length(ts_unemployment)), forecast_arima$lower[,2]),
Upper = c(rep(NA, length(ts_unemployment)), forecast_arima$upper[,2])
)
ggplot(df_arima, aes(x = Rok)) +
geom_line(aes(y = Actual), color = "pink") +
geom_line(aes(y = Predicted), color = "purple", linetype = "dashed") +
geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "purple", alpha = 0.2) +
labs(title = "ARIMA Forecast – unemployment (log)", y = "log-unemployment", x = "year") +
theme_minimal()## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_line()`).
# DIAGNOSTYKA MODEL_ARIMA
# Reszty modelu ARIMA
res_arima <- residuals(model_arima)
# Ljung-Box test (biały szum)
Box.test(res_arima, lag = 10, type = "Ljung-Box")##
## Box-Ljung test
##
## data: res_arima
## X-squared = 10.577, df = 10, p-value = 0.3914# Jarque-Bera Test (normalność reszt)
jarque.bera.test(res_arima)##
## Jarque Bera Test
##
## data: res_arima
## X-squared = 1.9192, df = 2, p-value = 0.3831# ARCH test (heteroskedastyczność)
ArchTest(res_arima)##
## ARCH LM-test; Null hypothesis: no ARCH effects
##
## data: res_arima
## Chi-squared = 2, df = 12, p-value = 0.9994# ACF i PACF reszt
par(mfrow = c(2, 2))
acf(res_arima, main = "ACF of ARIMA Residuals")
pacf(res_arima, main = "PACF of ARIMA Residuals")
hist(res_arima, main = "Histogram of Residuals", xlab = "Residuals", col = "lightblue", breaks = 10)
qqnorm(res_arima); qqline(res_arima, col = "red")
par(mfrow = c(1, 1))
# MIARY DOKŁADNOŚCI – MODEL_ARIMA
# ARIMA fitted values
fitted_arima <- fitted(model_arima)
actual_arima <- ts_unemployment
me_arima <- mean(fitted_arima - actual_arima)
rmse_arima <- sqrt(mean((fitted_arima - actual_arima)^2))
mae_arima <- mean(abs(fitted_arima - actual_arima))
mape_arima <- mean(abs((fitted_arima - actual_arima) / actual_arima)) * 100
mase_arima <- mean(abs(fitted_arima - actual_arima)) / mean(abs(diff(actual_arima)))
rmspe_arima <- sqrt(mean(((fitted_arima - actual_arima) / actual_arima)^2)) * 100
cat("\n== ARIMA Model Accuracy ==\n")##
## == ARIMA Model Accuracy ==cat("ME:", me_arima, "\nRMSE:", rmse_arima, "\nMAE:", mae_arima, "\nMAPE:", mape_arima,
"\nMASE:", mase_arima, "\nRMSPE:", rmspe_arima, "\n")## ME: -0.000184717
## RMSE: 0.1056651
## MAE: 0.08471971
## MAPE: 4.142828
## MASE: 0.7716067
## RMSPE: 5.306457Forecasting unemployment 2024-2030 based on the model
prognozytrend <- read_excel("C:/Users/zuzal/OneDrive/Desktop/prognozytrend.xlsx")
prognozy_ARIMA<- read_excel("C:/Users/zuzal/OneDrive/Desktop/prognozy_ARIMA.xlsx")
lata_prognozy <- 2024:2030
# Trend
# Dopasowane wartości (fitted) + prognoza
fitted_trend <- exp(fitted(model_log))
forecast_trend <- exp(predict(model_log, newdata = prognozytrend))
# Dane łączone: fitted + forecast
dane_trend <- data.frame(
rok = c(df$rok, lata_prognozy),
unemployment = c(fitted_trend, forecast_trend)
)
# ARIMA
prognoza_ue_zarimy <- predict(model_log, newdata = prognozy_ARIMA)
expv_arima <- exp(prognoza_ue_zarimy)
dane_arima <- data.frame(rok = Rok, unemployment = expv_arima)
# Do wykresu trendowego
dane_laczne_trend <- df %>%
select(rok, unemployment) %>%
mutate(typ = "Actual") %>%
bind_rows(dane_trend %>% mutate(typ = "Forecast (Trend)"))
# Do wykresu ARIMA
dane_laczne_arima <- df %>%
select(rok, unemployment) %>%
mutate(typ = "Actual") %>%
bind_rows(dane_arima %>% mutate(typ = "Forecast (ARIMA)"))
library(ggplot2)
# Wykres – TREND
ggplot(dane_laczne_trend, aes(x = rok, y = unemployment, color = typ)) +
geom_line(size = 1.2) +
labs(title = "Unemployment: trend forecast vs actual data",
x = "year", y = "unemployment") +
scale_color_manual(values = c("Actual" = "pink", "Forecast (Trend)" = "purple")) +
theme_minimal()
prognozytrend <- read_excel("C:/Users/zuzal/OneDrive/Desktop/prognozytrend.xlsx")
prognoza_ue_ztrendu <- predict(model_log, newdata=prognozytrend)
cat("Forecast result based on the model:", prognoza_ue_ztrendu,"\n")## Forecast result based on the model: 1.432686 1.345493 1.258265 1.171074 1.083874 0.996675 0.9094793expv_trend <- exp(prognoza_ue_ztrendu)
cat("Forecast result based on the model after de-logarithmization:", expv_trend)## Forecast result based on the model after de-logarithmization: 4.18994 3.84008 3.519309 3.225454 2.956108 2.709259 2.483029# Wykres – ARIMA
ggplot(dane_laczne_arima, aes(x = rok, y = unemployment, color = typ)) +
geom_line(size = 1.2) +
labs(title = "Unemployment: ARIMA forecast vs actual data",
x = "year", y = "unemployment") +
scale_color_manual(values = c("Actual" = "pink", "Forecast (ARIMA)" = "purple")) +
theme_minimal()
prognozyarima <- read_excel("C:/Users/zuzal/OneDrive/Desktop/prognozyarima.xlsx")
prognoza_ue_zarimy <- predict(model_log, newdata=prognozyarima)
expv_arima <- exp(prognoza_ue_zarimy)
cat("Forecast result based on the model:", prognoza_ue_zarimy,"\n")## Forecast result based on the model: 1.67946 1.759686 1.839917 1.920149 2.000408 2.080639 2.160871cat("Forecast result based on the model after de-logarithmization:", expv_arima)## Forecast result based on the model after de-logarithmization: 5.362661 5.810613 6.296016 6.821975 7.392074 8.009587 8.678695# RMSPE function
oblicz_rmspe <- function(rzeczywiste, dopasowane) {
sqrt(mean(((rzeczywiste - dopasowane) / rzeczywiste)^2)) * 100
}
# Dane rzeczywiste i dopasowane
rzeczywiste_2010_2023 <- df$unemployment
dopasowane_trend <- fitted_trend
dopasowane_arima <- fitted_arima
# RMSPE dla trendu
rmspe_trend <- oblicz_rmspe(rzeczywiste_2010_2023, dopasowane_trend)
# RMSPE dla ARIMA
rmspe_arima <- oblicz_rmspe(rzeczywiste_2010_2023, dopasowane_arima)
# Wyniki
cat("RMSPE for the Trend Model (2010–2023):", round(rmspe_trend, 2), "%\n")## RMSPE for the Trend Model (2010–2023): 4.52 %cat("RMSPE for the ARIMA Model (2010–2023):", round(rmspe_arima, 2), "%\n")## RMSPE for the ARIMA Model (2010–2023): 74.16 %SUMMARY
The forecasting of the unemployment rate in Poland for the period 2024-2030 initially required the selection of variables that could potentially influence the development of this indicator. Ultimately, five variables were identified.
In the process of model construction, all selected variables proved to be statistically significant. Subsequently, statistical tests were conducted to assess the adequacy and validity of the developed model.
The next phase of the study involved forecasting the explanatory variables for the years 2024-2030. Each variable was projected using two approaches: the linear trend method and the ARIMA (AutoRegressive Integrated Moving Average) method. Additionally, we tested other forecasting techniques, including the weighted moving average and Holt’s method; however, these models did not yield satisfactory results and were therefore not used in the final analysis.
In the case of the linear trend method, the variables were forecasted correctly and produced consistent results. However, when using the ARIMA model, the forecasts for expenditure of voivodships budgets and inflation did not perform well and deviated significantly from expected trends.
Also the linear trend method and ARIMA were applied again to generate the final unemployment forecasts. The values of the explanatory variables forecasted using both the linear trend method and ARIMA were then used to project the values of the unemployment rate using Linear Regression Model with logarithmized variables.
Finally, the unemployment rate values forecasted using both methods were compared with the actual unemployment data from the years 2010–2023, allowing us to evaluate the models’ performance and the reliability of the forecasts.
The forecast based on Linear Regression Model with logarithmized variables forcasted by linear trend method turned out to be the most accurate one. It’s RMSPE stands at 4.52 % which indicates that, on average, the model’s predictions deviate from the actual observed values by approximately 4.52%. This level of error suggests a very good fit of the model to the historical data, meaning that the model provides accurate and reliable forecasts for this variable. It means that the unemployment rate values will be approximately 4.18994 3.84008 3.519309 3.225454 2.956108 2.709259 2.483029 in the years 2024-2030.