v1-double-edge.R

#clean 
rm(list = ls())
library(sfcr)
library(tidyverse)
library(networkD3)
library(dplyr)


##########################################
#Equations of the model
##########################################


open_eqs <- sfcr_set(
  #Output
  Y_N ~ C_N + Ic_N + Id_N + G_N + X_N - IM_N,
  Y_S ~ C_S + Ic_S + Id_S + G_S + X_S - IM_S,
  
  
  #Industrial Output for emissions
  YKc_N ~ qc_N * Kc_N,
  YKd_N ~ qd_N * Kd_N,
  YK_N ~ YKc_N + YKd_N,
  
  YKc_S ~ qc_S * Kc_S,
  YKd_S ~ qd_S * Kd_S,
  YK_S ~ YKc_S + YKd_S,
  
  #productivity growth
  qc_N ~ qc_N[-1] + qc1_N + qc2_N * ggc_N + qc3_N * etac_N,
  qd_N ~ qd_N[-1] + qd1_N + qd2_N * ggd_N + qd3_N * etad_N,
  qc_S ~ qc_S[-1] + qc1_S + qc2_S * ggc_S + qc3_S * etac_S,
  qd_S ~ qd_N[-1] + qd1_S + qd2_S * ggd_S + qd3_S * etad_S,
  

  #accumulation rate
  gg_N ~ ((Kd_N + Kc_N)-(Kd_N[-1] + Kc_N[-1]) / Kd_N + Kc_N),
  gg_S ~ ((Kd_S + Kc_S)-(Kd_S[-1] + Kc_S[-1]) / Kd_S + Kc_S),

  ggc_N ~ (Ic_N / Kc_N[-1]),
  ggc_S ~ (Ic_S / Kc_S[-1]),
  ggd_N ~ (Id_N / Kd_N[-1]),
  ggd_S ~ (Id_S / Kd_S[-1]),
  
  
  #share in of capital type
  etac_N ~ (Kc_N/(Kd_N + Kc_N)),
  etad_N ~ (Kd_N/(Kd_N + Kc_N)),
  etac_S ~ (Kc_S/(Kd_S + Kc_S)),
  etad_S ~ (Kd_S/(Kd_S + Kc_S)),
  
  ####################
  #Households
  ####################
  #household income
  # consist of wages + profits + bank business profits
  YH_N ~ W_N + Pi_N + PiB_N,
  YH_S ~ W_S + Pi_S + PiB_S,
  # wages as a share of output
  W_N  ~ w_N * Y_N,
  W_S  ~ w_S * Y_S,
  #disposable income
  YD_N ~ YH_N - TX_N + r_N[-1] * Bh_N[-1],
  YD_S ~ YH_S - TX_S + r_S[-1] * Bh_S[-1],
 
  #wealth of HH - equation for behaviour
  V_N ~ V_N[-1] + ( YD_N - C_N ),
  V_S ~ V_S[-1] + ( YD_S - C_S ),
  
  #wealth of HH (balancing item)
  Vv_N ~ Hh_N + Bh_N + E_N,
  Vv_S ~ Hh_S + Bh_S + E_S,
  
  #household investment / equity demand
  Ih_N ~ lambda3_N * (V_N - Bh_N),
  Ih_S ~ lambda3_S * (V_S - Bh_S),
  
  #consumption
  C_N ~ alpha1_N * YD_N + alpha2_N * V_N[-1],
  C_S ~ alpha1_S * YD_S + alpha2_S * V_S[-1],
  #money held
  Hh_N ~ V_N - Bh_N - Ih_N,
  Hh_S ~ V_S - Bh_S - Ih_S,
  #portfolio choice
  Bh_N ~ V_N * ( lambda0_N + lambda1_N * r_N - lambda2_N * ( YD_N/V_N ) ),
  Bh_S ~ V_S * ( lambda0_S + lambda1_S * r_S - lambda2_S * ( YD_S/V_S ) ),
  
  
  ####################
  #Banks
  ####################
  #business profits of banks
  PiB_N ~ rd_N * Ld_N[-1] + rc_N * Lc_N[-1],
  PiB_S ~ rd_S * Ld_S[-1] + rc_S * Lc_S[-1],
  
  #wealth of banks
  Vcb_N ~ L_N,
  Vcb_S ~ L_S,
  

  
  ####################
  #Trade
  ####################
  
  #imports
  IM_N ~ mure_N * Y_N,
  IM_S ~ mu_S * Y_S,
  #responsive imports
  mure_N ~ mu_N + mu2_N * (Ic_N/Id_N),
  
  #change of imports
  deltamure_N ~ mure_N-mure_N[-1],
  
  #exports
  X_N ~ IM_S / xr,
  X_S ~ IM_N * xr,
  
 
  
  
  ####################
  #Firms
  ####################
  
  #firms profits
  Pi_N ~ C_N + Ic_N + Id_N + G_N - W_N + X_N - IM_N - rd_N * Ld_N[-1] - rc_N * Lc_N[-1],
  Pi_S ~ C_S + Ic_S + Id_S + G_S - W_S + X_S - IM_S - rd_S * Ld_S[-1] - rc_S * Lc_S[-1],
  

  # Investment
  Ic_N ~ Iac_N - I2c_N* ((rc_N*Xi_N)/rd_N) + I3c_N * (qc_N[-1]/qd_N[-1]),
  Id_N ~ Iad_N - I2d_N* (rd_N/(rc_N*Xi_N)) + I3d_N * (qd_N[-1]/qc_N[-1]),
  Ic_S ~ Iac_S - I2c_S* (rc_S/rd_S) + I3c_S * (qc_S[-1]/qd_S[-1]),
  Id_S ~ Iad_S - I2d_S* (rd_S/rc_S) + I3d_S * (qd_S[-1]/qc_S[-1]) + I4d_S * CRIT,
  
  #raw material demand
  CRIT ~ kappa * (Kc_N-Kc_N[-1]),

  #equity
  E_N ~ E_N[-1] + Ih_N,
  E_S ~ E_S[-1] + Ih_S,
  #kapital
  K_N ~ Kd_N + Kc_N,
  K_S ~ Kd_S + Kc_S,

  Kd_N ~ Kd_N[-1] + Id_N[-1] - deltad_N,
  Kc_N ~ Kc_N[-1] + Ic_N[-1] - deltac_N,
  Kd_S ~ Kd_S[-1] + Id_S[-1] - deltad_S,
  Kc_S ~ Kc_S[-1] + Ic_S[-1] - deltac_S,
  
  Kdg_N ~ d(Kd_N),
  Kcg_N ~ d(Kc_N),
  Kdg_S ~ d(Kd_S),
  Kcg_S ~ d(Kc_S),
  
  
  #Wealth of Firms
  Vf_N ~ Vf_N[-1] + Kd_N + Kc_N - Ld_N - Lc_N - E_N,
  Vf_S ~ Vf_S[-1] + Kd_S + Kc_S - Ld_S - Lc_S - E_S,
  
  
  #Loans
  L_N ~ Ld_N + Lc_N,
  L_S ~ Ld_S + Lc_S,
  
  
  Ld_N ~ Ld_N[-1] + Id_N - lambda4_N * Ih_N,
  Lc_N ~ Lc_N[-1] + Ic_N - (1-lambda4_N) * Ih_N,
  
  Ld_S ~ Ld_S[-1] + Id_S - lambda4_S * Ih_S,
  Lc_S ~ Lc_S[-1] + Ic_S - (1-lambda4_S) * Ih_S,
  
  ################
  #Government
  ################
  
  #taxation
  TX_N ~ theta_N * ( YH_N + r_N[-1] * Bh_N[-1] ),
  TX_S ~ theta_S * ( YH_S + r_S[-1] * Bh_S[-1] ),
  
  #Budget constraint of government
  Bs_N ~ Bs_N[-1] + ( G_N + r_N[-1] * Bs_N[-1] ) - ( TX_N + r_N[-1] * Bcb_N[-1] ),
  Bs_S ~ Bs_S[-1] + ( G_S + r_S[-1] * Bs_S[-1] ) - ( TX_S + r_S[-1] * Bcb_S[-1] ),
  #CB clearing of outstanding bills
  Bcb_N ~ Bs_N - Bh_N,
  Bcb_S ~ Bs_S - Bh_S,
  #balance sheet contraint of CB (via gold reserves)
  or_N ~ or_N[-1] + (( Hs_N - Hs_N[-1] - ( Bcb_N - Bcb_N[-1] ) )/pg_N),
  or_S ~ or_S[-1] + (( Hs_S - Hs_S[-1] - ( Bcb_S - Bcb_S[-1] ) )/pg_S),
  #supply of cash money
  Hs_N ~ Hh_N,
  Hs_S ~ Hh_S,
  #price of gold
  pg_S ~ pg_N * xr,
  #change of gold reserves
  deltaor_S ~ or_S - or_S[-1],
  deltaor_N ~ - (or_N - or_N[-1]),
  or ~ or_N + or_S,
  
  #interest rates
  #depend in interest rate on bonds could be endogenised with equity debt ratio
  rc_N ~ r_N,# * Xi_N,
  rd_N ~ r_N,
  rc_S ~ r_S,
  rd_S ~ r_S,
  
  
  #)
  
  #--- EMISSIONS BLOCK ---
  
  # Industrial CO₂ emissions (in GtCO2), output expressed in Northern currency
  EMIS_IN_N ~ CI_N  * YK_N,
  EMIS_IN_S ~ CI_S  * (YK_S  / xr),
  
  
  # Total (production‐based) emissions
  EMIS_TOTAL ~ EMIS_IN_N + EMIS_IN_S,
  
  # Carbon intensities (GtCO2 per unit of output in own currency)

  CI_N ~ CI_max - ((CI_max - CI_min) / (1 + ci_1 * exp(-ci_2 * (Kc_N[-1] / Kd_N[-1])))),
  CI_S ~ CI_max - ((CI_max - CI_min) / (1 + ci_1 * exp(-ci_2 * (Kc_S[-1] / Kd_S[-1])))),
  
  # Emissions growth rates and CI ratios, unchanged
  EMIS_growth_N ~ ((EMIS_IN_N - EMIS_IN_N[-1]) / EMIS_IN_N[-1]) * 100,
  EMIS_growth_S ~ ((EMIS_IN_S - EMIS_IN_S[-1]) / EMIS_IN_S[-1]) * 100,
  EMIS_growth_Total ~ ((EMIS_TOTAL - EMIS_TOTAL[-1]) / EMIS_TOTAL[-1]) * 100,
  CI_ratio_N   ~ Kc_N / Kd_N,
  CI_ratio_S   ~ Kc_S / Kd_S,
  
)

##########################################
#Parameters - all values fit for stability
##########################################

open_ext <- sfcr_set(
  xr ~ 4.447, #calibrated Exchange rate
  w_N ~ 0.6,
  w_S ~ 0.6,
  pg_N ~ 1,
  r_N ~ 0.025,
  r_S ~ 0.025,
  G_S ~ 20,
  G_N ~ 20,
  #Im and Export parameters
  mu_N ~ 0.15,
  mu_S ~ 0.35,
  mu2_N ~ 0.05,
  #mu2_S ~ 0.15, just for consistency
  #Consumption parameters
  alpha1_N ~ 0.7,
  alpha1_S ~ 0.7,
  alpha2_N ~ 0.3,
  alpha2_S ~ 0.3,
  #various size and allocation parameters
  lambda0_N ~ 0.67,
  lambda0_S ~ 0.67,
  lambda1_N ~ 0.05,
  lambda1_S ~ 0.05,
  lambda2_N ~ 0.01,
  lambda2_S ~ 0.01,
  lambda3_N ~ 0.006,
  lambda3_S ~ 0.006,
  lambda4_N ~ 0.5,
  lambda4_S ~ 0.5,
  theta_N ~ 0.2,
  theta_S ~ 0.2,
  #parameters of investment
  Iac_N ~ 0.25,
  Iad_N ~ 0.25,
  Iac_S ~ 0.25,
  Iad_S ~ 0.25,
  I1c_N ~ 0.1,
  I2c_N ~ 0.2,
  I3c_N ~ 0.022,
  I1d_N ~ 0.1,
  I2d_N ~ 0.2,
  I3d_N ~ 0.022,
  I1c_S ~ 0.1,
  I2c_S ~ 0.2,
  I3c_S ~ 0.022,
  I1d_S ~ 0.1,
  I2d_S ~ 0.2,
  I3d_S ~ 0.022,
  I4d_S ~ 0.5,
  kappa ~ 0.8,
  
  ###productivity - values taken from Quirino de Souza and Costa da Silva
  #kaldor-mirelees effect / constant
  qc1_N ~ 0.006,
  qd1_N ~ 0.006,
  qc1_S ~ 0.006,
  qd1_S ~ 0.006,
  #verdoorn coefficient 
  qc2_N ~ 0.05,
  qd2_N ~ 0.05,
  qc2_S ~ 0.05,
  qd2_S ~ 0.05,
  #structuralist share of capital in all kapital
  qc3_N ~ 0.02,
  qd3_N ~ 0.02,
  qc3_S ~ 0.02,
  qd3_S ~ 0.02,
  
  #kapital depreciation
  deltac_N  ~ 0.05,
  deltad_N  ~ 0.05,
  deltac_S  ~ 0.05,
  deltad_S  ~ 0.05,
  #green bond effect - this lowers the interest rate on green loans
  Xi_N ~ 1,
  Xi_S ~ 1,
 
  #=== EMISSIONS PARAMETERS from Dafermos & Nikolaidi ===
  CI_max ~ 0.003,      # Maximum carbon intensity (GtCO2/USD trillion) - Values are from Define Simple divided by 2 and divided by trillions (we simulate full values)
  CI_min ~ 0.00025,     # Minimum carbon intensity (GtCO2/USD trillion)
  ci_1 ~ 2.451037,   # Parameter linking green/conventional capital ratio to CI
  ci_2 ~ 3.579244    # Parameter linking green/conventional capital ratio to CI
  
)
########################
#setting initial values
# all other values start at 1
########################
open_initial <- sfcr_set(
  
  EMIS_IN_N ~ 0.0183,  # Initial emissions for North (GtCO2)
  EMIS_IN_S ~ 0.0183,   # Initial emissions for South (GtCO2)
  
  
)


# Simulate without hidden equations
open <- sfcr_baseline(
  equations = open_eqs,
  external = open_ext,
  initial = open_initial,
  periods = 100,
  hidden = c("deltaor_S" = "deltaor_N"
  ), .hidden_tol = 1)


#plotting to long format and adding CB and deficit developments

to_long <- function(model) {
  model %>%
    mutate(TB_N = X_N - IM_N,
           TB_S = X_S - IM_S,
           GB_N = TX_N - (G_N + dplyr::lag(r_N) * dplyr::lag(Bh_N)),
           GB_S = TX_S - (G_S + dplyr::lag(r_S) * dplyr::lag(Bh_S)),
           deltaV_N = V_N - lag(V_N),
           deltaV_S = V_S - lag(V_S),
           deltaBcb_N = Bcb_N - lag(Bcb_N),
           deltaBcb_S = Bcb_S - lag(Bcb_S),
           deltaHs_N = Hs_N - lag(Hs_N),
           deltaHs_S = Hs_S - lag(Hs_S)) %>%
    pivot_longer(cols = -period) %>%
    mutate(country = case_when(
      str_detect(name, '_N') ~ "North",
      str_detect(name, "_S") ~ "South",
      T ~ NA_character_
    ))
}

open_long <- to_long(open)


##########################
#Green bond shock
#########################

#Shocking the policy variable

shock1 <- sfcr_shock(
  variables = sfcr_set(Xi_N ~ 0.8), #Baseline setting Xi_N ~ 1 - for shock its 0.8
  start = 15,
  end = 30 #original value 30
)


open_1 <- sfcr_scenario(open, scenario = shock1, periods = 100)
open_1l <- to_long(open_1)


#Industry Formation in growth rates
open_1l %>%
  filter(name %in% c("Kcg_N", "Kdg_N", "Kcg_S", "Kdg_S")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  facet_wrap(~country, axes = "all", axis.labels = "all_y") +
  labs(title = "Industry Formation in growth rates")

#Twin Deficits
open_1l %>%
  filter(name %in% c("TB_S", "TB_N", "GB_S", "GB_N")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  facet_wrap(~country, axes = "all", axis.labels = "all_y") +
  labs(title = "Twin Deficits")

#Evolution of Exports and Imports
open_1l %>%
  filter(name %in% c("X_S", "X_N", "IM_N", "IM_S")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  #facet_wrap(~country, axes = "all", axis.labels = "all_y") +
  labs(title = "Change of Ex / Imports")

#Productivity growth
open_1l %>%
  filter(name %in% c("qc_N", "qd_N", "qc_S", "qd_S")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  facet_wrap(~country, axes = "all", axis.labels = "all_y") +
  labs(title = "Change of Productivity")


# Plot carbon emissions by region and globally
open_1l %>%
  filter(name %in% c("EMIS_IN_N", "EMIS_IN_S", "EMIS_TOTAL")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  labs(title = "Carbon Emissions Evolution",
       subtitle = "Industrial CO2 Emissions by Region and Global Total")

#Plot Carbon Intensity evolution in the North and South
open_1l %>%
  filter(name %in% c("CI_N", "CI_S")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  labs(title = "Carbon Emissions Intensity",
       subtitle = "Industrial CO2 Intensity by Region")

#Emissions growth rates in the North and South 

open_1l %>%
  filter(name %in% c("EMIS_growth_N", "EMIS_growth_S", "EMIS_growth_Total")) %>%
  ggplot(aes(x = period, y = value, color = name)) +
  geom_line(aes(linetype = name)) +
  labs(title = "Carbon Emissions Growth rates")



#in 65x15 a pdf export it is readable
sfcr_dag_cycles_plot(open_eqs)

library(ggraph)
library(tidygraph)

cycles_graph <- sfcr_dag_cycles(open_eqs)

ggraph(cycles_graph, layout = "kk") +
  geom_edge_link(arrow = arrow(type = "closed", length = unit(1, "mm")),
                 end_cap = circle(1, 'mm'), color = "grey") +
  geom_node_point(size = 1, color = "lightpink") +
  geom_node_text(aes(label = name), repel = TRUE, size = 3, color = "black" ) +
  theme_void()