#Toy model: moneta e circuito reddito-spesa
rm(list=ls())

#Parametri
theta=0.2 #0.3                    #Aliquota fiscale 
c_1=0.6                           #Propensione al consumo
T=25                              #Numero di periodi della simulaizone
lambda=0                          #% dei risparmi detenuti sotto forma di titoli pubblici
t_shock=15                        #Periodo in cui viene implementato lo shock (variazione spesa pubblica o aliquota fiscale o propensione al consumo) 
p=2
N=300
#Dichiarazione variabili
Y=matrix(data=0, nrow=T)          #Reddito lordo = PIL
C=matrix(data=0, nrow=T)          #Consumo delle famiglie
G=matrix(data=0, nrow=T)          #Spesa pubblica
Tax=matrix(data=0, nrow=T)        #Entrate fiscali
M=matrix(data=0, nrow=T)          #Depositi delle famiglie
B_h=matrix(data=0, nrow=T)        #Titoli pubblici delle famiglie
B_CB=matrix(data=0, nrow=T)       #Titoli pubblici detenuti dalla banca centrale
S=matrix(data=0, nrow=T)          #Risparmi delle famiglie
B=matrix(data=0, nrow=T)          #Debito pubblico
Deficit=matrix(data=0, nrow=T)    #Deficit pubblico
deltaS=matrix(data=0, nrow=T)     #Variazione dei risparmi in ogni periodo
Deficit_PIL=matrix(data=0,nrow=T) #Rapporto deficit pubblico-PIL 
I=matrix(data=0,nrow=T)           #Investimenti autonomi
Prestiti=matrix(data=0,nrow=T)      #Prestiti
L=matrix(data=0,nrow=T)           #Stock di debito delle imprese
Prestiti_PIL=matrix(data=0,nrow=T) 
Employment=matrix(data=0,nrow=T) 
deltaB_CB=matrix(data=0,nrow=T) 
deltaM=matrix(data=0,nrow=T) 
deltaBh=matrix(data=0,nrow=T) 
YD=matrix(data=0,nrow=T) 
for(t in 2:T){
  G[t]=50
  I[t]=0 #30
  if(t>=t_shock){                 # lo shock viene inserito a partire del periodo t_shock
    #G[t]=80
    #c_1=0.6
    #theta=0.3
    #I[t]=50
    }
  
  C[t]=Y[t-1]*(1-theta)*c_1 
  Y[t]=C[t]+G[t]+I[t]
  Employment[t]=Y[t]/p
  Tax[t]=Y[t]*theta
  Deficit[t]=G[t]-Tax[t]
  Deficit_PIL[t]=Deficit[t]/Y[t]
  B[t]=max(0,B[t-1]+G[t]-Tax[t])
  S[t]=S[t-1]+Y[t-1]*(1-theta)-C[t]
  deltaS[t]=S[t]-S[t-1]
  Prestiti[t]=I[t]
  L[t]=L[t-1]+Prestiti[t]
  Prestiti_PIL[t]=Prestiti[t]/Y[t]
  deltaBh[t]=min(deltaS[t]*lambda,max(Deficit[t],0))
  B_h[t]=B_h[t-1]+deltaBh[t]
  deltaM[t]=deltaS[t]- deltaBh[t]
  M[t]=M[t-1]+deltaM[t]
  deltaB_CB[t]=max(0,Deficit[t]-deltaBh[t])
  B_CB[t]=B_CB[t-1]+deltaB_CB[t]
  YD[t]=Y[t]*(1-theta)
}

a=3
b=T

#Layout grafico
m=rbind(c(1,2,3),c(4,5,6),c(7,8,9),c(10,11,12))
layout(m)
par(tcl=-0.1,mgp=c(2,0.4,0))
par(mar = c(1, 1.5, 1.5, 1), oma = c(0, 0, 1, 0))  # Aggiungi margine esterno superiore per i titoli

#GDP
minimo=min(Y[a:b])
massimo=max(Y[a:b])*1.1
plot(type='l',Y[a:b],lwd=3, col=4,ylim=range(minimo,massimo),ylab="",xlab="",main="GDP")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = Y[t_shock-1], y1 = Y[t_shock-1], col = 1, lwd = 1, lty = 2)

#Deficit Pubblico
minimo=min(Deficit[a:b])*0.6
massimo=max(Deficit[a:b])*1.2
plot(type='l',Deficit[a:b],lwd=3, col=4,ylim=range(minimo,massimo),ylab="",xlab="",main="Deficit Pubblico")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = Deficit[t_shock-1], y1 = Deficit[t_shock-1], col = 1, lwd = 1, lty = 2)


#Rapporto deficit/PIL
minimo=min(Deficit_PIL[a:b])*0.6
massimo=max(Deficit_PIL[a:b])*1
plot(type='l',Deficit_PIL[a:b],lwd=3, col=4,ylim=range(minimo,massimo),ylab="",xlab="",main="Deficit Pubblico/PIL")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = Deficit_PIL[t_shock-1], y1 = Deficit_PIL[t_shock-1], col = 1, lwd = 1, lty = 2)


#Variazione dei risparmi
minimo=min(deltaS[a:b])
massimo=max(deltaS[a:b])*1.1
plot(type='l',deltaS[a:b],lwd=3, col=3,ylim=range(minimo,massimo),ylab="",xlab="",main="Delta risparmi privati")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = deltaS[t_shock-1], y1 = deltaS[t_shock-1], col = 1, lwd = 1, lty = 2)

#Prestiti dele imprese
minimo=min(Prestiti[a:b])*0.9
massimo=max(Prestiti[a:b])*1.1
plot(type='l',Prestiti[a:b],lwd=3, col=3,ylim=range(minimo,massimo),ylab="",xlab="",main="Prestiti imprese")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = Prestiti[t_shock-1], y1 = Prestiti[t_shock-1], col = 1, lwd = 1, lty = 2)

#Prestiti privati/PIL
minimo=min(Prestiti_PIL[a:b])*0.8
massimo=max(Prestiti_PIL[a:b])*1.1
plot(type='l',Prestiti_PIL[a:b],lwd=3, col=3,ylim=range(minimo,massimo),ylab="",xlab="",main="DeltaL/PIL")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = Prestiti_PIL[t_shock-1], y1 = Prestiti_PIL[t_shock-1], col = 1, lwd = 1, lty = 2)


#Titoli pubblici acquistati dalla Banca Centrale
minimo=min(deltaB_CB[a:b])*0.8
massimo=max(deltaB_CB[a:b])*1.1
plot(type='l',deltaB_CB[a:b],lwd=3, col=6,ylim=range(minimo,massimo),ylab="",xlab="",main="Delta CB bonds")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = deltaB_CB[t_shock-1], y1 = deltaB_CB[t_shock-1], col = 1, lwd = 1, lty = 2)

#Variazione dei depositi delle famiglie
minimo=min(deltaM[a:b])*0.8
massimo=max(deltaM[a:b])*1.1
plot(type='l',deltaM[a:b],lwd=3, col=6,ylim=range(minimo,massimo),ylab="",xlab="",main="Delta Deposits")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = deltaM[t_shock-1], y1 = deltaM[t_shock-1], col = 1, lwd = 1, lty = 2)

#Titoli pubblici acquistati dalle famiglie
minimo=min(deltaBh[a:b])*0.8
massimo=max(deltaBh[a:b])*1.1
plot(type='l',deltaBh[a:b],lwd=3, col=6,ylim=range(minimo,massimo),ylab="",xlab="",main="Delta HH bonds")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)
segments(x0 = t_shock-3, x1 = b, y0 = deltaBh[t_shock-1], y1 = deltaBh[t_shock-1], col = 1, lwd = 1, lty = 2)



#Debito Pubblico
minimo=min(B[a:b])
massimo=max(B[a:b])*1.1
plot(type='l',B[a:b],lwd=3, col=7,ylim=range(minimo,massimo),ylab="",xlab="",main="Debito pubblico")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)


#Stock di risparmi delle famiglie
minimo=min(S[a:b])
massimo=max(S[a:b])*1.1
plot(type='l',S[a:b],lwd=3, col=7,ylim=range(minimo,massimo),ylab="",xlab="",main="Risparmi privati")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)

#Stock di depositi
minimo=min(M[a:b])
massimo=max(M[a:b])*1.1
plot(type='l',M[a:b],lwd=3, col=7,ylim=range(minimo,massimo),ylab="",xlab="",main="Depositi")
polygon(c(t_shock-3, rep(b, 2), t_shock-3), c(minimo, minimo, massimo, massimo), col= rgb(0,0,1.0,alpha=0.02), border = NA)