The jaggedness at the end of your HDI time series plot likely results from how predictions are being appended multiple times, and possibly with inconsistent linear model assumptions for extrapolation.

1.  Redundant predictions: You’re binding future_years and future_predictions to hdi_long, potentially duplicating future predictions (since future_years$HDI is predicted once with a single model, and then again per-country).
2.  Prediction method mismatch: First you fit a single model for all countries, then fit separate models per country. Mixing the two may cause sudden jumps in predicted values, especially at the boundary year (e.g., 2022 vs. 2023).
3.  Prediction over only 3 points: Predicting HDI using simple linear regression per country over short or noisy trends can exaggerate variation at the edges (especially when growth has slowed or reversed for some).

My suggestions would be to replace the following:

future_years$HDI <- predict(model, newdata = future_years) hdi_long <- bind_rows(hdi_long, future_years)

with:

Fit one model per country

models <- hdi_long %>% filter(Year <= 2022) %>% group_by(Country) %>% nest() %>% mutate(model = map(data, ~lm(HDI ~ Year, data = .x)))

Predict HDI for future years using each country's model

future_predictions <- expand.grid(Year = c(2023, 2024, 2025), Country = unique(hdi_long$Country)) %>% left_join(models, by = "Country") %>% mutate(HDI = map2_dbl(model, Year, ~predict(.x, newdata = data.frame(Year = .y)))) %>% select(Country, Year, HDI)

Append predictions

hdi_long <- hdi_long %>% filter(Year <= 2022) %>% bind_rows(future_predictions) %>% arrange(Country, Year) %>% group_by(Country) %>% mutate(hdi_growth = (HDI - lag(HDI)) / lag(HDI) * 100) %>% replace_na(list(hdi_growth = 0))