etfvsfp.py

# -*- coding: utf-8 -*-
"""ETFvsFP.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1H6B4APvArNIE83ZYt8LQ-rqitHcE89zO
"""

# Commented out IPython magic to ensure Python compatibility.
# tfr https://www.finanzaonline.com/forum/etf-fondi-e-gestioni-e-investment-certificates/1979835-chi-vi-da-un-rendimento-annuale-pari-alla-deducibilita-di-un-fondo-pensione-post56606644.html#post56606644
# tfr https://www.finanzaonline.com/forum
# simulazioni https://www.finanzaonline.com/forum/etf-fondi-e-gestioni-e-investment-certificates/1979835-chi-vi-da-un-rendimento-annuale-pari-alla-deducibilita-di-un-fondo-pensione-post56488416.html#post56488416
import pandas as pd
import numpy as np
import bottleneck as bn
import matplotlib.pyplot as plt
ARR_COMM = 2


# Dati lavoratore
timeHorizon = 40
dt = pd.date_range(start='2021', periods=timeHorizon, freq='y')
years = dt.year
RAL = 33000.0
premioStraordinario = 0.0
str_regione = 'Lombardia'
str_comune = 'Brescia'
# Dati TFR
TFR_percentuale = 1/13.5
TFR_impostaSostitutiva = 0.17 # https://www.agenziaentrate.gov.it/portale/web/guest/schede/pagamenti/f24-imposta-sostitutiva-tfr/cosa-f24-imposta-sostitutiva-tfr-imprese
FOI = np.full(shape=(timeHorizon,), fill_value=0.03, dtype=float)
# Dato INPS
INPS_sogliaPrimaFascia = 48279
INPS_sogliaMassima = 105014
INPS_soglie = [INPS_sogliaPrimaFascia, INPS_sogliaMassima]
INPS_percentuali = [0.0949, 0.01]
# Dati IRPEF
IRPEF_scaglioni = np.array([15000, 28000, 50000], dtype=np.float64)
IRPEF_aliquote = np.array([0.23, 0.25, 0.35, 0.43])
IRPEF_aliquoteAddizionaleRegionale = np.array([0.0123, 0.0158, 0.0172, 0.0173])
#IRPEF_scaglioniAddizionaleComunale = np.array([23000], dtype=np.float64) # Milano
IRPEF_scaglioniAddizionaleComunale = np.array([13000], dtype=np.float64) # Brescia
IRPEF_aliquoteAddizionaleComunale = np.array([0, 0.008]) # uguale sia a Milano che Brescia
# Dati FP
FP_taxMax = 0.15
FP_taxMin = 0.09
FP_taxDecr = 0.003 # decremento tassa annuale
FP_taxStartDecr = 15 # dopo quanti anni inizia decrescita tassazione
FPquotaDeducibile = 5165
#FPpercentualeDatore = 0.0155 # Fonte
#FPpercentualeLavoratore = 0.0055 # Fonte
FPpercentualeDatore = 0.022 # Cometa 2% ma 2.2% per nuove adesioni https://www.contrattometalmeccanici.it/art-15-previdenza-complementare
FPpercentualeLavoratore = 0.012 # Cometa
#FPpercentualeDatore = 0.01
#FPpercentualeLavoratore = 0.01
#FP_str = 'Espero Bilanciato'
FP_str = 'Cometa Bilanciato'
FP_costoGestione = np.full(shape=(timeHorizon,), fill_value=22, dtype=float)
#FP_returnRate = np.full(shape=(timeHorizon,), fill_value=0.0389, dtype=float) # Espero
#FP_returnRate = np.full(shape=(timeHorizon,), fill_value=0.04, dtype=float) # generico
FP_returnRate = np.full(shape=(timeHorizon,), fill_value=0.0498, dtype=float) # Cometa, comparto Crescita, media ultimi 10 anni
FPA_returnRate = np.full(shape=(timeHorizon,), fill_value=0.05, dtype=float) # Allianz Insieme Linea Azionaria (Azionario)
#FPA_returnRate = np.full(shape=(timeHorizon,), fill_value=0.0725, dtype=float) # Allianz Insieme Linea Azionaria (Azionario)
#FPA_returnRate = np.full(shape=(timeHorizon,), fill_value=0.0273, dtype=float) # Amundi SecondaPensione Espansione ESG (Azionario)
'''
GARC_W0GE = np.array([(825.67/100)**(1/37) - 1]) # 85' - 22'
GARC_G0Q0 = np.array([(700.424/100)**(1/37) - 1]) # x' - 22'
GARC_ER00 = np.array([(300.348/100)**(1/37) - 1]) # x' - 22'
GARC_W1SGITRE = np.array([0.0940]) # 12' - 22'
GARC_DJSWICLN = np.array([0.0715]) # 12' - 22'
GARC_FP_FonteDinamico = 0.35*GARC_W1SGITRE + 0.25*GARC_DJSWICLN + 0.25*GARC_W0GE + 0.05*GARC_G0Q0 + 0.10*GARC_ER00
'''

# Dati ETF
ETF_Bollo_Tax = np.array([0.002])
ETF_CG_TAX = np.array([0.26])
ETF_returnRate = np.full(shape=(timeHorizon,), fill_value=0.0876, dtype=float) # ACWI '87 - 22

# Dataframe
anni = np.arange(start=1, stop=timeHorizon+1, dtype=np.uint)
RALs = np.full(shape=(timeHorizon,), fill_value=RAL, dtype=float)
#RALs = np.repeat(np.array([6, 9, 10, 11, 11]), 3)*5000
isCambioAzienda = np.full_like(RALs, False, dtype=bool)
#isCambioAzienda[[10, 20, 30]] = True
#isCambioAzienda[[3, 6, 9, 12, 15, 20]] = True
RALs = np.linspace(RAL, RAL+20000, timeHorizon, dtype=float)

df = pd.DataFrame(RALs, columns=['RAL'], index=years)

df.head()

# Sezione INPS
df['premio_straordinario'] = premioStraordinario
df['imponibile_INPS'] = df['RAL'] + df['premio_straordinario']

# Implementazione logica if-elif-else vettoriale
'''
def calcoloContributiINPS(imponibileINPS, sogliaPrimaFascia, sogliaMassimale, IVS_PERC, IVS_AGG_PERC):
  if imponibileINPS <= sogliaPrimaFascia:
      IVS = IVS_PERC * imponibileINPS
      IVS_agg = 0.0
  elif imponibileINPS > sogliaPrimaFascia and imponibileINPS <= sogliaMassimale:
      IVS = IVS_PERC * imponibileINPS
      IVS_agg = (imponibileINPS - sogliaPrimaFascia)*IVS_AGG_PERC
  else: # imponibileINPS > sogliaMassimale
      IVS = IVS_PERC * sogliaMassimale
      IVS_agg = (sogliaMassimale - sogliaPrimaFascia)*IVS_AGG_PERC
  return [IVS, IVS_agg, IVS+IVS_agg]
'''
# 2021
def calcola_IVS(imponibileINPS):
  # if-elif
  condlist = [imponibileINPS <= INPS_soglie[0],
              (imponibileINPS > INPS_soglie[0]) & (imponibileINPS <= INPS_soglie[1])]
  choicelist = [INPS_percentuali[0]*imponibileINPS,
                INPS_percentuali[0]*imponibileINPS]
  # else
  default = INPS_percentuali[0]*INPS_soglie[1]
  return np.select(condlist, choicelist, default)

# 2021
def calcola_IVS_aggiuntiva(imponibileINPS):
  # if-elif
  condlist = [imponibileINPS <= INPS_soglie[0],
              (imponibileINPS > INPS_soglie[0]) & (imponibileINPS <= INPS_soglie[1])]
  choicelist = [0.,
                (imponibileINPS-INPS_soglie[0]) * INPS_percentuali[1]]
  # else
  default = (INPS_soglie[1]-INPS_soglie[0]) * INPS_percentuali[1]
  return np.select(condlist, choicelist, default)

df['IVS'] = calcola_IVS(df['imponibile_INPS'].to_numpy())
df['IVS'] = df['IVS'].round()

df['IVS_aggiuntiva'] = calcola_IVS_aggiuntiva(df['imponibile_INPS'].to_numpy())
df['IVS_aggiuntiva'] = df['IVS_aggiuntiva'].round()

df['contributo_INPS'] = df['IVS'] + df['IVS_aggiuntiva']

def calcolaImponibileIRPEF(imponibileINPS, contributoINPS, FPquotaVolontaria):
  return imponibileINPS - contributoINPS - FPquotaVolontaria

# vectorized function:
# logica if elif else tramite combinazione lineare
def calcolaScaglioniProgressivi(x, scaglioni, aliquote):
  idx = np.searchsorted(scaglioni, x)
  #print(idx[0])
  scaglioniPrependZero = np.zeros((aliquote.size,))
  scaglioniPrependZero[1:] = scaglioni
  #print(scaglioniPrependZero)
  # logica:
  # first = np.diff(scaglioniPrependZero[:idx[0]+1]) @ aliquote[:idx[0]]
  # last = (x-scaglioniPrependZero[idx[0]]) * aliquote[idx[0]]
  # print(first)
  # print(last)
  # return first + last
  a = np.empty((idx[0]+1,))
  a[:-1] = np.diff(scaglioniPrependZero[:idx[0]+1])
  a[-1] = x-scaglioniPrependZero[idx[0]]
  b = aliquote[:idx[0]+1]
  return a @ b

'''
def calcolaImposta(imponibileIRPEF, IRPEF_scaglioni, IRPEF_aliquote):
  if imponibileIRPEF <= IRPEF_scaglioni[0]:
    imposta = 0     + IRPEF_aliquote[0] * (imponibileIRPEF-0);
  elif imponibileIRPEF > IRPEF_scaglioni[0] and imponibileIRPEF <= IRPEF_scaglioni[1]:
    imposta = 3450  + IRPEF_aliquote[1] * (imponibileIRPEF-IRPEF_scaglioni[0])
  elif imponibileIRPEF > IRPEF_scaglioni[1] and imponibileIRPEF <= IRPEF_scaglioni[2]:
    imposta = 6700  + IRPEF_aliquote[2] * (imponibileIRPEF-IRPEF_scaglioni[1])
  else:
    imposta = 14400 + IRPEF_aliquote[3] * (imponibileIRPEF-IRPEF_scaglioni[2])
  return round(imposta, ARR_COMM)
'''
def calcolaImpostaIRPEF(imponibileIRPEF, IRPEF_scaglioni, IRPEF_aliquote):
  return np.around(np.apply_along_axis(calcolaScaglioniProgressivi, 0, np.atleast_2d(imponibileIRPEF), IRPEF_scaglioni, IRPEF_aliquote), ARR_COMM)

# https://www.regione.lombardia.it/wps/portal/istituzionale/HP/DettaglioRedazionale/servizi-e-informazioni/cittadini/tributi-e-canoni/addizionale-irpef
# Lombardia
'''
def calcolaAddizionaleRegionale(imponibileIRPEF, IRPEF_scaglioni, IRPEF_aliquoteAddizionaleRegionale):
  if imponibileIRPEF < IRPEF_scaglioni[0]:
    addizionaleRegionale = IRPEF_aliquoteAddizionaleRegionale[0]*imponibileIRPEF
  elif imponibileIRPEF > IRPEF_scaglioni[0] and imponibileIRPEF <= IRPEF_scaglioni[1]:
    addizionaleRegionale = 184.5 + IRPEF_aliquoteAddizionaleRegionale[1]*(imponibileIRPEF-IRPEF_scaglioni[0])
  elif imponibileIRPEF > IRPEF_scaglioni[1] and imponibileIRPEF <= IRPEF_scaglioni[2]:
    addizionaleRegionale = 389.9 + IRPEF_aliquoteAddizionaleRegionale[2]*(imponibileIRPEF-IRPEF_scaglioni[1])
  else:
    addizionaleRegionale = 768.3 + IRPEF_aliquoteAddizionaleRegionale[3]*(imponibileIRPEF-IRPEF_scaglioni[2])
  addizionaleRegionale = round(addizionaleRegionale, ARR_COMM)
'''
def calcolaAddizionaleRegionale(imponibileIRPEF, IRPEF_scaglioni, IRPEF_aliquoteAddizionaleRegionale):
  return np.around(np.apply_along_axis(calcolaScaglioniProgressivi, 0, np.atleast_2d(imponibileIRPEF), IRPEF_scaglioni, IRPEF_aliquoteAddizionaleRegionale), ARR_COMM)

# https://www.comune.milano.it/aree-tematiche/tributi/addizionale-comunale-irpef
# Milano
'''
def calcolaAddizionaleComunale(self):
  if self.imponibile <= 23000:
    self.addizionaleComunale = 0
  else:
    self.addizionaleComunale = 0.008*self.imponibile 
  self.addizionaleComunale = round(self.addizionaleComunale, ARR_COMM) 
'''
def calcolaAddizionaleComunale(imponibileIRPEF, IRPEF_scaglioniAddizionaleComunale, IRPEF_aliquoteAddizionaleComunale):
  return np.around(np.piecewise(imponibileIRPEF,
                      [imponibileIRPEF <= IRPEF_scaglioniAddizionaleComunale[0]],
                      [lambda x: IRPEF_aliquoteAddizionaleComunale[0]*x,
                       lambda x: IRPEF_aliquoteAddizionaleComunale[1]*x]),
                   ARR_COMM)

def calcolaAliquotaIRPEF(imponibileIRPEF, impostaIRPEF):
  return impostaIRPEF/imponibileIRPEF

def calcolaAliquotaMediaMobileIRPEF(aliquotaIRPEF, window):
  # TODO sistemare media mobile nei primi cinque anni
  return bn.move_mean(aliquotaIRPEF, window, min_count=1)

#@detrazioneProQuota(365)
def calcolaDetrazioneIRPEF(imponibileIRPEF_DetrLavDip, IRPEF_scaglioni):
  return np.around(np.piecewise(imponibileIRPEF_DetrLavDip,
                      [imponibileIRPEF_DetrLavDip <= IRPEF_scaglioni[0],
                       (imponibileIRPEF_DetrLavDip > IRPEF_scaglioni[0]) & (imponibileIRPEF_DetrLavDip <= IRPEF_scaglioni[1]),
                       (imponibileIRPEF_DetrLavDip > IRPEF_scaglioni[1]) & (imponibileIRPEF_DetrLavDip <= IRPEF_scaglioni[2])],
                      [lambda x: 1880,
                       lambda x: 1910 + 1190 * ((IRPEF_scaglioni[1]-x)/(IRPEF_scaglioni[1]-IRPEF_scaglioni[0])),
                       lambda x:        1910 * ((IRPEF_scaglioni[2]-x)/(IRPEF_scaglioni[2]-IRPEF_scaglioni[1])),
                       lambda x: 0]),
                   ARR_COMM)
  
def calcolaDetrazioneExtraIRPEF(imponibileIRPEF_DetrLavDip):
  return ((imponibileIRPEF_DetrLavDip >= 25000) & (imponibileIRPEF_DetrLavDip <= 35000)) * 65

def checkNoTaxArea(detrazioneTotaleIRPEF, impostaTotaleIRPEF):
  # https://www.finanzaonline.com/forum/etf-fondi-e-gestioni-e-investment-certificates/1979835-chi-vi-da-un-rendimento-annuale-pari-alla-deducibilita-di-un-fondo-pensione-post56664207.html#post56664207
  # TODO: controllare condizione no tax area:
  # Metodo 1: self.imponibile <= 8174
  # Metodo 2: self.imposta <= 1880 
  return detrazioneTotaleIRPEF >= impostaTotaleIRPEF

def calcolaCapienzaDisponibile(impostaIRPEF, addizionaleRegionale, addizionaleComunale):
  return impostaIRPEF + addizionaleRegionale + addizionaleComunale
  
def calcolaCapienzaRichiesta(detrazioneIRPEF, bonusIRPEF):
  return detrazioneIRPEF + bonusIRPEF

def calcolaRedditoNetto(self):
  if self.isNoTaxArea:
    self.detrazioni = self.imposta
    self.redditoNetto = self.imponibile
  else:
    self.calcolaBonus()
    self.calcolaAddizionaleRegionale()
    self.calcolaAddizionaleComunale()
    self.calcolaCapienzaDisponibile()
    self.calcolaCapienzaRichiesta()
    if self.capienzaDisponibile < self.capienzaRichiesta:
      # possibile assorbire i bonus/derazioni solo se si può scontare dalle imposte (irpef e addizionali)
      # altrimenti reddito netto = reddito lordo (imponibile irpef)
      self.redditoNetto = self.imponibile
    else:
      self.redditoNetto = self.imponibile + self.capienzaRichiesta - self.capienzaDisponibile
    self.redditoNetto = round(self.redditoNetto, ARR_COMM)
    
# Bonus ex-Renzi
def calcolaBonus(imponibileIRPEF, IRPEF_scaglioni):
  return (imponibileIRPEF <= IRPEF_scaglioni[0]) * 1200

# Logica
def updateCalcola(imponibileINPS, FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico):
  FPquotaVolontaria = FPquotaVolontariaBustaPaga + FPquotaVolontariaBonifico
  imponibileIRPEF = calcolaImponibileIRPEF(imponibileINPS, df['contributo_INPS'].to_numpy(), FPquotaVolontaria)
  imponibileIRPEF_DetrLavDip = calcolaImponibileIRPEF(imponibileINPS, df['contributo_INPS'].to_numpy(), FPquotaVolontariaBustaPaga)
  impostaIRPEF = calcolaImpostaIRPEF(imponibileIRPEF, IRPEF_scaglioni, IRPEF_aliquote)
  detrazioneIRPEF = calcolaDetrazioneIRPEF(imponibileIRPEF_DetrLavDip, IRPEF_scaglioni)
  detrazioneExtraIRPEF = calcolaDetrazioneExtraIRPEF(imponibileIRPEF_DetrLavDip)
  detrazioneTotaleIRPEF = detrazioneIRPEF + detrazioneExtraIRPEF
  bonusIRPEF = calcolaBonus(imponibileIRPEF, IRPEF_scaglioni)
  isNoTaxArea = checkNoTaxArea(detrazioneTotaleIRPEF, impostaIRPEF)
  addizionaleRegionale = np.zeros_like(imponibileIRPEF)
  addizionaleComunale = np.zeros_like(imponibileIRPEF)
  impostaTotaleIRPEF = np.zeros_like(imponibileIRPEF)
  addizionaleRegionale[~isNoTaxArea] = calcolaAddizionaleRegionale(imponibileIRPEF[~isNoTaxArea], IRPEF_scaglioni, IRPEF_aliquoteAddizionaleRegionale)
  addizionaleComunale[~isNoTaxArea] = calcolaAddizionaleComunale(imponibileIRPEF[~isNoTaxArea], IRPEF_scaglioniAddizionaleComunale, IRPEF_aliquoteAddizionaleComunale)
  impostaTotaleIRPEF = impostaIRPEF + addizionaleRegionale + addizionaleComunale
  IRPEF_netta = np.clip(impostaTotaleIRPEF - detrazioneTotaleIRPEF, 0., None)
  redditoNetto = imponibileIRPEF - IRPEF_netta + bonusIRPEF
  aliquotaIRPEF = calcolaAliquotaIRPEF(imponibileIRPEF, impostaIRPEF)
  aliquotaMAIRPEF = calcolaAliquotaMediaMobileIRPEF(aliquotaIRPEF, window=5)

  cols = ['quota volontaria', 'quota volontaria (busta paga)', 'quota volontaria (bonifico)',
          'reddito complessivo', 'Imponibile DetrLavDip', 'aliquota media', 'aliquota media mobile',
          'IRPEF', 'addizionale regionale', 'addizionale comunale', 'detrazione', 'detrazione extra', 'bonus Renzi',
          'imposta totale', 'detrazione totale', 'IRPEF netta', 'reddito netto']

  dataIRPEF = np.column_stack([FPquotaVolontaria, FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico,
                              imponibileIRPEF, imponibileIRPEF_DetrLavDip, aliquotaIRPEF, aliquotaMAIRPEF,
                              impostaIRPEF, addizionaleRegionale, addizionaleComunale, detrazioneIRPEF, detrazioneExtraIRPEF, bonusIRPEF,
                              impostaTotaleIRPEF, detrazioneTotaleIRPEF, IRPEF_netta, redditoNetto])
  return pd.DataFrame(dataIRPEF, columns=cols)

imponibileTFR = RALs

# Contributi previdenza integrativa a carico datore
# https://www.finanzaonline.com/forum/etf-fondi-e-gestioni-e-investment-certificates/1979835-chi-vi-da-un-rendimento-annuale-pari-alla-deducibilita-di-un-fondo-pensione-post56834329.html#post56834329
# quota datoriale e quota aderante sono calcolate sull'imponibile TFR (o talvolta sul minimo minimo retributivo come Cometa)
minimoRetributivoMensile = 2003.99 # livello B1
minimoRetributivoAnnuale = minimoRetributivoMensile * 13 # CCNL Metalmeccanico ha 13 mensilità

#FPquotaDatoriale = np.around(FPpercentualeDatore * imponibileTFR, ARR_COMM) # Fonte
FPquotaDatoriale = np.around(FPpercentualeDatore * minimoRetributivoAnnuale, ARR_COMM) # Cometa (calcolata su minimi contrattuali)

FPquotaVolontariaMassima = FPquotaDeducibile - FPquotaDatoriale

# Contributi previdenza integrativa a carico lavoratore – Deducibili
FPquotaVolontariaMinima = np.around(FPpercentualeLavoratore * imponibileTFR, ARR_COMM)

# Caso A: Investo in un ETF
FPquotaVolontariaBustaPaga = np.zeros_like(RALs)
FPquotaVolontariaBonifico = np.zeros_like(RALs)
dfIRPEF_A = updateCalcola(df['imponibile_INPS'].to_numpy(), FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico) 
dfIRPEF_A[:5]

# Caso B: Investo in un FP con il contributo datoriale e saturando la deducibilità (FPN o FPA convenzionato)
FPquotaVolontariaBustaPaga = np.full_like(RALs, FPquotaDeducibile - FPquotaDatoriale)
FPquotaVolontariaBonifico = np.zeros_like(RALs)
dfIRPEF_B = updateCalcola(df['imponibile_INPS'].to_numpy(), FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico) 
dfIRPEF_B[:5]

# Caso C: Investo in un FP (FPN o FPA convenzionato) con il contributo minimo per avere il contributo datoriale e il resto in ETF 
FPquotaVolontariaBustaPaga = np.full_like(RALs, FPquotaVolontariaMinima)
FPquotaVolontariaBonifico = np.zeros_like(RALs)
dfIRPEF_C = updateCalcola(df['imponibile_INPS'].to_numpy(), FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico) 
dfIRPEF_C[:5]

# Caso D: Investo in un FP (FPN o FPA convenzionato) con il contributo minimo per avere il contributo datoriale e il resto in FPA
# considero solo FPA perché posso sfruttare il caso C per vedere l'altro contributo... (TODO FIX: non coincede con caso sotto)
FPquotaVolontariaBustaPaga = np.zeros_like(RALs)
FPquotaVolontariaBonifico = np.full_like(RALs, FPquotaDeducibile - (FPquotaVolontariaMinima + FPquotaDatoriale))
dfIRPEF_D = updateCalcola(df['imponibile_INPS'].to_numpy(), FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico) 
dfIRPEF_D[:5]

# Caso D2: Investo in un FP (FPN o FPA convenzionato) con il contributo minimo per avere il contributo datoriale e il resto in FPA
# (TODO FIX: non coincede con caso sopra)
FPquotaVolontariaBustaPaga = np.full_like(RALs, FPquotaVolontariaMinima)
FPquotaVolontariaBonifico = np.full_like(RALs, FPquotaDeducibile - (FPquotaVolontariaMinima + FPquotaDatoriale))
dfIRPEF_D2 = updateCalcola(df['imponibile_INPS'].to_numpy(), FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico) 
dfIRPEF_D2[:5]

# Caso F: Investo in un FPA
FPquotaVolontariaBustaPaga = np.zeros_like(RALs)
FPquotaVolontariaBonifico = np.full_like(RALs, FPquotaDeducibile)
dfIRPEF_F = updateCalcola(df['imponibile_INPS'].to_numpy(), FPquotaVolontariaBustaPaga, FPquotaVolontariaBonifico) 
dfIRPEF_F[:5]

Contributo_A_ETF = np.around(dfIRPEF_A.loc[:, 'reddito netto'].to_numpy() - dfIRPEF_B.loc[:, 'reddito netto'].to_numpy(), ARR_COMM)
Contributo_A_ETF_effettivo = Contributo_A_ETF

Contributo_B_FP = np.full((timeHorizon,), FPquotaDeducibile, dtype=float)
Contributo_B_FP_effettivo = Contributo_A_ETF

Contributo_C_FP = np.full((timeHorizon,), FPquotaDatoriale + FPquotaVolontariaMinima, dtype=float)
Contributo_C_FP_effettivo = np.around(dfIRPEF_A.loc[:, 'reddito netto'].to_numpy() - dfIRPEF_C.loc[:, 'reddito netto'].to_numpy(), ARR_COMM)
Contributo_C_ETF =  np.around(Contributo_A_ETF_effettivo - Contributo_C_FP_effettivo, ARR_COMM)
Contributo_C_ETF_effettivo =  Contributo_C_ETF

Contributo_D_FP = Contributo_C_FP
Contributo_D_FP_effettivo = Contributo_C_FP_effettivo
Contributo_D_FPA = np.full((timeHorizon,), FPquotaDeducibile - FPquotaVolontariaMinima - FPquotaDatoriale, dtype=float)
Contributo_D_FPA_effettivo = np.around(dfIRPEF_A.loc[:, 'reddito netto'].to_numpy() - dfIRPEF_D.loc[:, 'reddito netto'].to_numpy(), ARR_COMM)

Contributo_E_ETF = Contributo_D_FPA_effettivo + Contributo_D_FP_effettivo
Contributo_E_ETF_effettivo = Contributo_E_ETF

Contributo_F_FPA = np.full((timeHorizon,), FPquotaDeducibile, dtype=float)
Contributo_G_ETF = np.around(dfIRPEF_A.loc[:, 'reddito netto'].to_numpy() - dfIRPEF_F.loc[:, 'reddito netto'].to_numpy(), ARR_COMM)

Contributo_A_ETF_effettivo[0]

Contributo_B_FP_effettivo[0]

Contributo_C_FP_effettivo[0]+Contributo_C_ETF_effettivo[0]

Contributo_E_ETF_effettivo[0]

Contributo_D_FP_effettivo[0]+Contributo_D_FPA_effettivo[0]

Contributo_D2_effettivo = np.around(dfIRPEF_A.loc[:, 'reddito netto'].to_numpy() - dfIRPEF_D2.loc[:, 'reddito netto'].to_numpy(), ARR_COMM)
Contributo_D2_effettivo[0]

# (A-D+A-C) - (A-B) = 2A-C-D-A+B = A+B-C-D
# A-D2 - (A-B) = B - D2
costoBonifico = Contributo_E_ETF_effettivo[0]-Contributo_A_ETF_effettivo[0]
costoBonifico

Contributo_D2_effettivo[0] - Contributo_B_FP_effettivo[0]

def costoCompraVenditaETF(controvalore):
  return np.around(np.clip(controvalore*.0019, 1.5, 5.), ARR_COMM)

def tassaCapitalFP(year):
  return np.clip(FP_taxMax - FP_taxDecr*(year - FP_taxStartDecr), FP_taxMin, FP_taxMax)

def simulazioneETF(A_gross, returnRate):
  finalValue = np.zeros((timeHorizon,))
  A_net = np.zeros((timeHorizon,))
  bollo = np.zeros((timeHorizon,))
  capital = np.zeros((timeHorizon,))
  finalValue_gross = np.zeros((timeHorizon,), dtype=float)
  FV = 0
  for tt in range(timeHorizon):
    A_net[tt] = A_gross[tt] - costoCompraVenditaETF(A_gross[tt])
    FV = FV + A_net[tt]
    finalValue_gross[tt] = FV
    bollo[tt] =  np.around(ETF_Bollo_Tax * FV, ARR_COMM)
    FV = FV * (1. + returnRate[tt])
    FV = np.around(FV, ARR_COMM)
    # il bollo non lo sottraggo al capitale investito perché
    # non vendo quote per pagarlo ma suppongo di versare liquidità in più
  capital = np.cumsum(A_gross)
  bolloSum = np.cumsum(bollo)
  capitalGain = finalValue_gross - capital
  capitalGainTaxes = np.around(capitalGain * ETF_CG_TAX, ARR_COMM)
  finalValue_net = finalValue_gross - capitalGainTaxes - bolloSum
  return np.column_stack([finalValue_net, finalValue_gross, capitalGainTaxes, capitalGain, bollo, bolloSum, capital])

FP_anniContribuzionePrecedente = 0.

def simulazioneFP(A_gross, returnRate, costoGestione, timeHorizon):
  A_net = np.zeros((timeHorizon,))
  capital = np.zeros((timeHorizon,))
  finalValue_gross = np.zeros((timeHorizon,))
  FV = 0
  for tt in range(timeHorizon):
    A_net[tt] = A_gross[tt] - costoGestione[tt]
    FV = FV + A_net[tt]
    finalValue_gross[tt] = FV
    FV = FV * (1. + returnRate[tt])
    FV = np.around(FV, ARR_COMM)
  capital = np.cumsum(A_gross)
  capitalTaxes = np.around(capital * tassaCapitalFP(np.arange(1, timeHorizon+1))+FP_anniContribuzionePrecedente, ARR_COMM)
  finalValue_net = finalValue_gross - capitalTaxes
  return np.column_stack([finalValue_net, finalValue_gross, capitalTaxes, capital])


def simulazioneFPalt(A_gross, returnRate, costoGestione, timeHorizon):
  #A_net = np.zeros((timeHorizon,))
  #capital = np.zeros((timeHorizon,))
  #finalValue_gross = np.zeros((timeHorizon,))
  #finalValue_net = np.zeros((timeHorizon,))
  #res_alt = np.column_stack([finalValue_net, finalValue_gross, capitalTaxes, capital])
  res_alt = np.zeros((timeHorizon, 4))
  idxs = np.flatnonzero(isCambioAzienda)
  idxs_diff = np.diff(idxs, prepend=0, append=timeHorizon)
  idx = 0
  for ii in range(idxs_diff.size):
    FP_anniContribuzionePrecedente = idx
    res = simulazioneFP(A_gross[idx:idx+idxs_diff[ii]], returnRate, costoGestione, idxs_diff[ii])
    if ii > 0:
      res_alt[idx:idx+idxs_diff[ii], :] = res + res_alt[idx-1, :]
    else:
      res_alt[idx:idx+idxs_diff[ii], :] = res

    if ii == idxs_diff.size-1:
      break
    else:
      idx = idxs[ii]
  return res_alt

'''
def renditaAnticipata(R, i, t):
  return R*(1+i)*((1+i)**t-1)/i

def renditaPosticipata(R, i, t):
  return R*((1+i)**t-1)/i

def calcoloTassaCapitaleAnticipata(b, R, i, t):
  return b*R*(1+i)*((((1+i)**(t+1)-1)/i)-(t+1))/i

def calcoloTassaCapitalePosticipata(b, R, i, t):
  return b*R*((((1+i)**(t+1)-1)/i)-(t+1))/i
'''

# Sezione TFR
def returnRateTFR(FOI, impostaSostitutiva):
  return (0.015 + 0.75 * FOI) * (1 - impostaSostitutiva)

def simulazioneAziendaTFR(A_gross, returnRate, timeHorizon, aliquotaMediaMobileCinqueAnni):
  finalValue = np.zeros((timeHorizon,))
  capital = np.zeros((timeHorizon,))
  finalValue_gross = np.zeros((timeHorizon,), dtype=float)
  FV = 0
  for tt in range(timeHorizon):
    FV = FV + A_gross[tt]
    finalValue_gross[tt] = FV
    FV = FV * (1. + returnRate[tt])
    FV = np.around(FV, ARR_COMM)
  capital = np.cumsum(A_gross)
  capitalTaxes = np.around(capital * aliquotaMediaMobileCinqueAnni, ARR_COMM)
  finalValue_net = finalValue_gross - capitalTaxes
  return np.column_stack([finalValue_net, finalValue_gross, capitalTaxes, capital])
  
def simulazioneAziendaTFRalt(A_gross, returnRate, timeHorizon, aliquotaMediaMobileCinqueAnni):
  #A_net = np.zeros((timeHorizon,))
  #capital = np.zeros((timeHorizon,))
  #finalValue_gross = np.zeros((timeHorizon,))
  #finalValue_net = np.zeros((timeHorizon,))
  #res_alt = np.column_stack([finalValue_net, finalValue_gross, capitalTaxes, capital])
  res_alt = np.zeros((timeHorizon, 4))
  idxs = np.flatnonzero(isCambioAzienda)
  idxs_diff = np.diff(idxs, prepend=0, append=timeHorizon)
  idx = 0
  for ii in range(idxs_diff.size):
    res_alt[idx:idx+idxs_diff[ii], :] = simulazioneAziendaTFR(A_gross[idx:idx+idxs_diff[ii]], returnRate, idxs_diff[ii], aliquotaMediaMobileCinqueAnni[idx:idx+idxs_diff[ii]])
    if ii == idxs_diff.size-1:
      break
    else:
      idx = idxs[ii]
  return res_alt

contributoTFR = np.around(imponibileTFR * TFR_percentuale, ARR_COMM)
# https://www.agenziaentrate.gov.it/portale/web/guest/schede/pagamenti/f24-imposta-sostitutiva-tfr/cosa-f24-imposta-sostitutiva-tfr-imprese
TFR_impostaSostitutiva = 0.17
TFR_returnRate = returnRateTFR(FOI, TFR_impostaSostitutiva)
aliquotaMediaMobileCinqueAnni = dfIRPEF_A.loc[:, 'aliquota media mobile'].to_numpy()


Montante_TFR_Azienda = simulazioneAziendaTFR(contributoTFR, TFR_returnRate, timeHorizon, aliquotaMediaMobileCinqueAnni)

Montante_TFR_AziendaAlt = simulazioneAziendaTFRalt(contributoTFR, TFR_returnRate, timeHorizon, aliquotaMediaMobileCinqueAnni)



Montante_TFR_AziendaAlt[np.flatnonzero(isCambioAzienda)-1, 0]

Montante_TFR_FP = simulazioneFP(contributoTFR, FP_returnRate, np.zeros_like(RALs), timeHorizon)

Montante_A_ETF = simulazioneETF(Contributo_A_ETF, ETF_returnRate)


Contributo_A2_ETF = Contributo_A_ETF
Contributo_A2_ETF[np.flatnonzero(isCambioAzienda)] += Montante_TFR_AziendaAlt[np.flatnonzero(isCambioAzienda)-1, 0]
Montante_A2_ETF = simulazioneETF(Contributo_A2_ETF, ETF_returnRate)

Montante_B_FP = simulazioneFP(Contributo_B_FP, FP_returnRate, FP_costoGestione, timeHorizon)

Montante_C_FP = simulazioneFP(Contributo_C_FP, FP_returnRate, FP_costoGestione, timeHorizon)

Montante_C_ETF = simulazioneETF(Contributo_C_ETF, ETF_returnRate)

Montante_D_FP = simulazioneFP(Contributo_D_FP, FP_returnRate, FP_costoGestione, timeHorizon)

Montante_D_FPA = simulazioneFP(Contributo_D_FPA, FPA_returnRate, FP_costoGestione, timeHorizon)

Montante_E_ETF = simulazioneETF(Contributo_E_ETF, ETF_returnRate)

Montante_F_FPA = simulazioneETF(Contributo_F_FPA, FPA_returnRate)

Montante_G_ETF = simulazioneETF(Contributo_G_ETF, ETF_returnRate)

# np.full(shape=(timeHorizon,), fill_value=0.0273, dtype=float) # Amundi SecondaPensione Espansione ESG (Azionario)
Montante_F2_FPA = simulazioneETF(Contributo_F_FPA, np.full(shape=(timeHorizon,), fill_value=0.0273, dtype=float))

ms = 2
x = np.arange(1, timeHorizon+1)
fig, ax = plt.subplots(figsize=(20,15))
plt.yscale('log')
plt.plot(x, Montante_B_FP[:, 0], marker='o', markersize=ms, label='Caso B: '+'FP='+str(Contributo_B_FP[0]))
plt.plot(x, Montante_A_ETF[:, 0], marker='o', markersize=ms, label='Caso A: '+'ETF='+str(Contributo_A_ETF[0]))
plt.plot(x, Montante_C_FP[:, 0]+Montante_C_ETF[:, 0], marker='o', markersize=ms, label='Caso C: '+'FP='+str(Contributo_C_FP[0])+' ,ETF='+str(Contributo_C_ETF[0]))
#plt.plot(x, Montante_D_FP[:, 0]+Montante_D_FPA[:, 0], marker='o', markersize=ms, label='Caso D: '+'FP='+str(Contributo_D_FP[0])+' ,FPA='+str(Contributo_D_FPA[0]))
#plt.plot(x, Montante_E_ETF[:, 0], marker='o', markersize=ms, label='Caso E: '+'ETF='+str(Contributo_E_ETF[0]))
plt.title(f'RAL: {RAL}')
plt.xlabel('Time [years]')
plt.ylabel('Final value [eur]')
plt.xticks(np.arange(1,timeHorizon+1,1))
plt.xlim(1, timeHorizon+1)
plt.grid(True, 'minor')
plt.legend();

ms = 2
x = np.arange(1, timeHorizon+1)
fig, ax1 = plt.subplots(figsize=(20,15))
plt.yscale('log')
idx = np.flatnonzero(isCambioAzienda)
for kk in idx:
  plt.axvline(kk, color='k', linestyle='--')
plt.plot(x, Montante_A_ETF[:, 0] + Montante_TFR_Azienda[:, 0], marker='o', markersize=ms, label='Caso A: '+'ETF='+str(Contributo_A_ETF[0])+', TFR(AZ)='+str(contributoTFR[0]))
plt.plot(x, Montante_B_FP[:, 0] + Montante_TFR_FP[:, 0], marker='o', markersize=ms, label='Caso B: '+'FP='+str(Contributo_B_FP[0])+' ('+str(Contributo_A_ETF[0])+'), TFR(FP)='+str(contributoTFR[0]))
plt.plot(x, Montante_C_FP[:, 0] + Montante_C_ETF[:, 0] + Montante_TFR_FP[:, 0], marker='o', markersize=ms, label='Caso C: '+'FP='+str(Contributo_C_FP[0])+' ('+str(Contributo_C_FP_effettivo[0])+'), ETF='+str(Contributo_C_ETF[0])+', TFR(FP)='+str(contributoTFR[0]))
#plt.plot(x, Montante_A2_ETF[:, 0] + Montante_TFR_AziendaAlt[:, 0], marker='o', markersize=ms, label='Caso A2: '+'ETF='+str(Contributo_A_ETF[0])+', TFR(AZ)='+str(contributoTFR[0]))
plt.title(f'({str_regione}/{str_comune}), rendimento ETF: {ETF_returnRate[0]}, rendimento FP: {FP_returnRate[0]}, FOI: {FOI[0]}')
plt.xlabel('Time [years]')
plt.ylabel('Final value [eur]')
plt.xticks(np.arange(1, timeHorizon+1))
plt.xlim(1, timeHorizon)
plt.ylim(2*10**3, 3*10**6)
plt.grid(True, 'minor')
plt.legend(loc='upper left', frameon=True, framealpha=1);

# 1. Choose your desired colormap
cmap = plt.get_cmap('Set2')
# 2. Segmenting the whole range (from 0 to 1) of the color map into multiple segments
slicedCM = cmap(np.linspace(0, 1, 8)) 

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
plt.plot(x, RALs, label='RAL', color=slicedCM[0])
#plt.plot(x, Montante_TFR_Azienda[:, 0], marker='o', markersize=ms, label='TFR Azienda', color=slicedCM[0])
plt.plot(x, Montante_TFR_AziendaAlt[:, 0], marker='o', markersize=ms, label='TFR Azienda',color=slicedCM[1])
plt.plot(x, Montante_TFR_FP[:, 0], marker='o', markersize=ms, label='TFR Fondo Pensione', color=slicedCM[2])

plt.tick_params(axis='y')
plt.ylim(0, 100000);
plt.legend(loc='upper center', frameon=True, framealpha=1);
#fig.tight_layout()  # otherwise the right y-label is slightly clipped

Montante_B_FP_alt = simulazioneFPalt(Contributo_B_FP, FP_returnRate, FP_costoGestione, timeHorizon)
Montante_B_FP_alt[-5:]

ms = 2
x = np.arange(1, timeHorizon+1)
fig, ax1 = plt.subplots(figsize=(20,15))
plt.yscale('log')
plt.plot(x, Montante_B_FP_alt[:, 0], marker='o', markersize=ms, label='Caso B: '+'FP='+str(Contributo_B_FP[0]))
plt.plot(x, Montante_B_FP[:, 0], marker='o', markersize=ms, label='Caso B: '+'FP='+str(Contributo_B_FP[0]))

plt.title(f'({str_regione}/{str_comune}), rendimento FP: {FP_returnRate[0]}')
plt.xlabel('Time [years]')
plt.ylabel('Final value [eur]')
idx = np.flatnonzero(isCambioAzienda)
for kk in idx:
  plt.axvline(kk, color='k')
plt.xticks(np.arange(1, timeHorizon+1))
plt.xlim(1, timeHorizon)
plt.ylim(2*10**3, 3*10**6)
plt.grid(True, 'minor')
plt.legend();
ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
plt.plot(x, RALs, label='RAL')
plt.tick_params(axis='y')
plt.ylim(0, 100000);
plt.legend();