'Optimal Debt Targets
'José Federico Geli - jfgeli@gmail.com
'José Marín Arcas - jose.marin@airef.es
'Ver 1.0
'21/Dec/2016

'********************************************PURPOSE*********************************************
'This program simulates default probabilities (i.e. probability of reaching a given debt target)
'by applying stochastic shocks to the variables.

call settings
call create_wf

call parameters
call seeder
call model_definition
call model_solution("_mod","d_r bal_r")

'********************************************SETTINGS*********************************************
subroutine settings 
mode quite 
logmode all
cd "E:\Google Drive\Temp\Optimal debt targets"

'solve options

!obs_n=200
!smpl_start = 2016
endsub

'********************************************WORKING ENVIRONMENT*********************************************
subroutine create_wf
!smpl_end = !smpl_start+!obs_n-1
'wfcreate(wf=ODT, page=Data) a {!smpl_start} {!smpl_end}  '1000 PERIODS
setmaxerrs 2
wfclose ODT
clearerrs
setmaxerrs 1
wfcreate(wf=ODT, page=Data) a 1900 {!smpl_end}  '1000 PERIODS
model _mod
'starting values 
smpl @all

	'Sample for estimation (Eviews format)
%esmpl="1900 2016"
' Sample for forecasting (Eviews format) 
%fsmpl="2017 2100"
' Variable for fan chart
%varname="d_r"
'Number of bands in fanchart starting from 90%. 
'For example, !numb=1 plots 90% band and the mode, 
'!numb=2 plots 90% band, 45% band and the mode 
!numb=4
endsub

'********************************************PREDEFINED PARAMETERS*********************************************
subroutine parameters
'Pre-defined parameters
!mu=0.12
!sigma=0.06
!r_sens = 0.05
!d_r_t=0.60
!d_r_market=0.90
!u=0.05
!v=0.5
!switch_t=1
!switch_g=0
!switch_i=0
!switch_fatigue=1


'Exogenous variables and parameters
!eta=0.005 'population growth rate
!epsilon=0.015 'labour growth rate
!alfa=0.7 'CD elasticity of consumption
!pi=0.02 'CPI evolution
!r=0.035 'debt interest rate
!d_start= 0.825
!r_elast = 0.05

'Shocks to implement
%ShockList = "eta epsilon alfa pi"
scalar eta_s=1' 0.02 'population growth rate
scalar epsilon_s=1 ' 0.01 'labour growth rate
scalar alfa_s=0 'CD elasticity of consumption
scalar pi_s=0 'CPI evolution

'Exogenous initial levels
!n=1 'number of people  
!k=1 'productivity of labour
!tau=0.385
!g_level=0.147010877361483
!i_level=0.18
!p=1 'price level
!shock_size=0.02

endsub

subroutine seeder
'********************************************SERIES SEED*********************************************


series shock = !shock_size*@runif(-1,1)

series u=!u
series v=!v
series eta=!eta+eta_s*shock
series epsilon=!epsilon+epsilon_s*shock
series alfa = !alfa+alfa_s*shock
series p = 1
smpl 2017 @last
series p=p(-1)*(1+!pi+pi_s*shock) 
smpl @all
series n=!n
n.label(d) Workers
series k=!k
k.label(d) Capital
series tau=!tau
series r = !r
series r_elast = !r_elast
series d_r_t = !d_r_t
series g = -1 + (1+eta)*(1+epsilon)*(1+0.02)
g.label(d) Nominal output growth

smpl 1900 2016
for %aInstrument {%ShockList}
	if {%aInstrument}_s=1 then
		equation {%aInstrument}_eq.ls {%aInstrument} c {%aInstrument}(-1)
	endif 
next 
smpl @all



smpl {!smpl_start} {!smpl_start}

series bal_pri_r_t = ((r-g)/(1+g))*d_r_t
bal_pri_r_t.label(d) Primary balance target
series g_level=!g_level
g_level.label(d) Productive public spending
series i_level=!i_level
i_level.label(d) Redistributive public spending
series  leis= (1-alfa)*(n+(g_level/(p*k))+(i_level/(p*k*(1-tau)))) 'leisure 
leis.label(d) Leisure
series cons= alfa*(1-tau)*(n*k+(g_level/p)+(i_level/(p*(1-tau)))) 'consumption
cons.label(d) Consumption
series y=k*(n-leis)+g_level/p 'real output level 
y.label(d) Real output
series y_level = y
y_level.label(d) Nominal output 

series i_r=i_level/y_level
i_r.label(d) Redistributive public spending/GDP ratio
series g_r=g_level/y_level
g_r.label(d) Productive public spending/GDP ratio
series d_level=!d_start-(!tau*y_level-g_level-i_level-r*!d_start)
d_level.label(d) Debt 
series d_r = d_level/y_level
d_r.label(d) Debt/GDP ratio
series  inp=r*!d_start 
inp.label(d) interest payments 

inp.label(d) interest payments 
series bal = tau*y_level-g_level-i_level-r*!d_start
bal.label(d) Budget balance
series bal_pri =tau*y_level-g_level-i_level 'primary balance 
bal_pri.label(d) Primary balance
series bal_r = bal/y_level
bal_r.label(d) Balance/GDP ratio
series  inp_r=r*!d_start/y_level
inp_r.label(r) interest payments/GDP ratio
series bal_pri_r = bal_pri/y_level
bal_pri_r.label(d) Primary balance/GDP ratio
series expen = inp+g_level+i_level
expen.label(d) Total expenditure
series expen_r = expen/y_level
expen_r.label(d) Total expenditure/GDP ratio
'series u_level = (cons^alfa)*(leis^(1-alfa)) 'utility level
'u_level.label(d) Utility
series y_disp = y_level*(1-tau)+i_level
y_disp.label(d) Disposable income

smpl @all
endsub

'********************************************ECONOMY MODEL*********************************************
subroutine model_definition
'levels requiring a starting value
_mod.append g=-1 + (1+eta)*(1+epsilon)*(1+@pc(p)/100) 'nominal output growth
_mod.append bal_pri_r_t = ((r-g)/(1+g))*d_r_t 'target primary balance ratio
_mod.append bal_pri_d = @recode((u(-1)*(d_level(-1)/y_level(-1)-d_r_t(-1))-v(-1)*(bal_pri(-1)/y_level(-1)-bal_pri_r_t(-1)))+bal_pri(-1)/y_level(-1)<0.03,(u(-1)*(d_level(-1)/y_level(-1)-d_r_t(-1))-v(-1)*(bal_pri(-1)/y_level(-1)-bal_pri_r_t(-1))),0.03-bal_pri(-1)/y_level(-1)) 'u,v rule, change in primary balance
_mod.append n=n(-1)*(1+eta) 'workers
_mod.append k=k(-1)*(1+epsilon)  'capital
_mod.append tau=tau(-1)+!switch_t*bal_pri_d
_mod.append g_level=(g_level(-1)/y_level(-1)-!switch_g*bal_pri_d)*(1+eta)*(1+epsilon)*(1+@pc(p)/100)*y_level(-1)
_mod.append i_level=(i_level(-1)/y_level(-1)-!switch_i*bal_pri_d)*(1+!eta)*(1+!epsilon)*(1+@pc(p)/100)*y_level(-1)


'endogenous levels
_mod.append bal=tau*y_level-g_level-i_level-r*d_level(-1) 'deficit
_mod.append inp=r*d_level(-1) 
_mod.append bal_pri=tau*y_level-g_level-i_level 'primary deficit 
_mod.append d_level = d_level(-1)-bal 'debt level
_mod.append leis= (1-alfa)*(n+(g_level/(p*k))+(i_level/(p*k*(1-tau)))) 'leisure 
_mod.append cons= alfa*(1-tau)*(n*k+(g_level/p)+(i_level/(p*(1-tau)))) 'consumption
_mod.append y=k*(n-leis)+g_level/p 'real output level 
_mod.append r=@recode(r(-1)+r_elast*(d_level(-1)/y_level(-1)-!d_r_market)>0,r_elast*(d_level(-1)/y_level(-1)-!d_r_market)+!r,!r)
'_mod.append r=!r
_mod.append y_level = y_level(-1)*(y/y(-1))*(p/p(-1)) 'nominal output level
_mod.append p=p(-1)*(1+!pi) 'price level
'_mod.append u_level = (cons^alfa)*(leis^(1-alfa)) 'utility level
_mod.append i_r=i_level/y_level
_mod.append g_r=g_level/y_level
_mod.append y_disp = y_level*(1-tau)+i_level
_mod.append d_r = d_level/y_level
_mod.append bal_r = bal/y_level
_mod.append bal_pri_r = bal_pri/y_level
_mod.append expen = inp+g_level+i_level
_mod.append expen_r = expen/y_level

for %aInstrument {%ShockList}
	if {%aInstrument}_s=1 then
		_mod.merge {%aInstrument}_eq 
	endif 
next 

endsub 
'_mod.solveopt(s=d,d=d,o=b)
'smpl 2016+1 @last
'_mod.solve 
'group _results bal_r_0 bal_pri_r_0 r_0 d_r_0
'group _output g_r inp_r g
'graph _results_graph.line  _results
'_results_graph.setelem(4) axis(r)
'smpl 2016 @last
'show _results_graph
'show bal_0 d_level_0 leis_0 cons_0 y_level_0 @pc(p_0) r_0
'call fan_charts


subroutine fan_charts
	
'	'Sample for estimation (Eviews format)
'%esmpl="1900 2016"
'
'' Sample for forecasting (Eviews format) 
'%fsmpl="2017 2050"

'Bootstrap errors? 1=yes, 0=no, use normal distribution
!boot=0

'Coefficient uncertainty? 1=yes, 0=no coefficient uncertanty
!coefu=0



'Do not change the 2 parameters below unless you are adapting the code
' for your own model
'Name of the model to be used for constructing a fan chart
%modname="_mod"

'Variable to forecast (has to be one of the model variables)
'%varname="d_r"
	
	
			smpl {%esmpl}
'=========== Model specification =========
'replace the two lines below with your model of inflation if you want
'equation eq1.ls {%varname} c {%varname}(-1) '{%varname}(-2)
'eq1.makemodel({%modname}) 
'=============== Main Fan Chart Program==============
%quant_array="90 "
!quant_step=90/!numb
!k=90-!quant_step
while !k>5 
%quant_array=%quant_array+@str(@round(!k))+" "
!k=!k-!quant_step
wend
%quant_array=%quant_array+"1"


delete(noerr) g2
group g2

for %quant {%quant_array}  
	!alfa={%quant} 'current confidence level
	smpl {%fsmpl}
	!alfa_ratio=!alfa/100 
	
	!num_simul=100

	if @abs(!boot-1)<0.0001 then 'boot strap errors
		if @abs(!coefu-1)<0.0001 then 'coefficient uncertainty
			{%modname}.stochastic(i=b, b=!alfa_ratio,c=t,r=!num_simul)
		else 'no coefficient uncertainty
			{%modname}.stochastic(i=b, b=!alfa_ratio,r=!num_simul)
		endif
	else 'normal errors
		if @abs(!coefu-1)<0.0001 then 'coefficient uncertainty
			{%modname}.stochastic(i=n, b=!alfa_ratio,c=t,r=!num_simul)
		else 'no coefficient uncertainty
			{%modname}.stochastic(i=n, b=!alfa_ratio,r=!num_simul)
		endif
	endif	
	

	{%modname}.solveopt(s=b,d=d)
	{%modname}.solve
	
	%first_elem_fsmpl=@wleft(%fsmpl,1)
	%second_elem_fsmpl=@wright(%fsmpl,1)
	
	smpl {%first_elem_fsmpl} {%first_elem_fsmpl}
	{%varname}_0h={%varname}
	{%varname}_0l={%varname}
	
	'delete(noerr) smpl_gr
	'sample smpl_gr {%first_elem_fsmpl}-6 {%second_elem_fsmpl}
	smpl {%first_elem_fsmpl}-6 {%second_elem_fsmpl}
	if !alfa<3 then
		%quant="0"
	endif
	series {%varname}_{%quant}h_={%varname}_0h
	series {%varname}_{%quant}l_={%varname}_0l
	g2.add {%varname}_{%quant}l_ {%varname}_{%quant}h_

next
series p1={%varname}
g2.add {%varname} p1
freeze(fanchart, mode=overwrite) g2.band(o="modern")

fanchart.legend -display
for !i=1 to !numb+1
	!r=@round(0.2*255)
	!b=0
	!g=@round((1-0.175*@sqrt(!i) )*255)
	fanchart.setelem(!i) fcolor(@rgb(!r,!g,!b))
next
!numline=!numb+2
fanchart.setelem(!numline) fcolor(blue) lwidth(10)


fanchart.draw(shade,b, gray) {%fsmpl}
show fanchart

smpl @all

endsub

subroutine model_solution(string %modname,string %trackedList)

'%trackedList = "d_r bal_r"
%solutionList = @wreplace(%trackedList,"*","*_0")

smpl 2017 2066
	
{%modname}.solveopt(s=b,d=d)
{%modname}.stochastic(i=n, r=100,f=0.85,p=stacked_results)
{%modname}.track {%trackedList}
{%modname}.solve
pageselect stacked_results
pageunstack(namepat=*?, page=unstacked_results) repid obsid @ {%solutionList}
	
	
endsub