CAL 2014:1 1
allocate  1000

*******************************************************************************
***                              Question                                   ***
*******************************************************************************

/*

On this code:
The following is a code to replicate Aizenman et al (2000) which investigates
a relationship between an income shock and an optimal level of the government
bond issues. They consider a simple two period setting where an income shock in
period 2 is uncertain between [-limit, limit].
As is well-known, if we have a negative income shock in period 1, the government
will issue their bonds for consumption smoothing. However, in this model,
the borrowing constraint is set depending on tax revenue in period 2 which in turn
depends on the expected income shocks in period 2. Therefore, they may unable to
accomplish their optimazation perfectly. In a polar case, their bond issue will
be fixed at a certain level.The following code is to depict the situation above.
When you run it, you can see the figure the bond issue is restricted at a certain
level.

My question:
The revenue in period 2 depends on tax efficiency parameter lambda as well.
When lambda is, for example, 1.1, this code exactly replicate Aizenman (2000).
However, when lambda = 1.55, we see weird horizontal part of the graph which is not
emerged in Aizenman (2000). Actually, when I changed the method from Simplex to,
for example, to generic, I have a different figure again. Thus, I expect the reason
why I got the different figure from Aizenman (2000) is my inferior technique for Rats.
So, could you tell me if you have any idea? I am looking forward to your assistance.

Reference:
Joshua Aizenman , Michael Gavin & Ricardo Hausmann (2000) Optimal tax
and debt policy with endogenously imperfect creditworthiness, The Journal of International Trade &
Economic Development, 9:4, 367-395, DOI: 10.1080/096381900750056830

To link to this article:
https://doi.org/10.1080/096381900750056830

*/

*******************************************************************************
***                        Key parameter (tax efficiency)                   ***
*******************************************************************************
compute lambda  = 1.55

*******************************************************************************
***                              Initial settings                           ***
*******************************************************************************
compute Y1     = 1.0
compute Y2     = 1.0
compute G      = 0.2
compute rho    = 0.1
compute alpha  = 0.5
compute xi_1   = 0.2
compute C1     = 0.8
compute GrossRate = 1.1

compute RiskAversion  = 1.2
compute dU1       = 1.0
compute E_dU2     = 1.0
compute shock1    = 0.0
compute s_limit   = 0.3
compute threshold = 0.0
compute Bond      = 0.0

dec vec shock2(100)
dec vec Y2_shocked(100)
dec vec xi_2(100)
dec vec C2(100)
dec vec dU2(100)

compute estimation_upper = 35

dec vec E_dU21(estimation_upper)
dec vec shock1_for_figure(estimation_upper)
dec vec Bond_for_figure(estimation_upper)


*******************************************************************************
***                Loop of various level of shocks in period 1              ***
*******************************************************************************
do j = 1, estimation_upper
compute shock1 = -float(j)/100.0

********************
*** Optimization ***
********************
nonlin Bond
find(METHOD=simplex,iterations=100) root dU1 - E_dU2

****************************************
*** Gross Interest Rates in Period 1 ***
****************************************
compute Grossrate = %if((1.0+rho)*Bond<=alpha*((Y2-s_limit)/(2.0*lambda)-G),1.0+rho, $
alpha/Bond*((Y2+s_limit)/(2.0*lambda)-G) - (alpha*s_limit)/(lambda*Bond)*sqrt(((2.0*lambda)/s_limit)*((Y2/(2.0*lambda))-G) $
- 2.0*lambda*Bond*(1.0+rho)/(alpha*s_limit)))

************************************
*** Marginal Utility in Period 1 ***
************************************
compute xi_1    = (G-Bond)/(Y1+shock1)
compute C1      = (1.0 - xi_1-(1.0/(2.0*lambda))*(1 - sqrt(1 - 2.0*lambda*xi_1))^2)*(Y1+shock1)
compute dU1     = C1^(-RiskAversion)*(1.0/sqrt(1.0-2.0*lambda*xi_1))

*********************************************
*** Expected Marginal Utility in Period 2 ***
*********************************************
compute threshold  = (2.0*lambda)*((Grossrate*Bond)/alpha + G)  - Y2
compute lowerb  = %max(-s_limit, threshold)
do i = 1, 100
compute shock2(i) = (s_limit - lowerb)*float(i)/100.0 + lowerb
compute Y2_shocked(i) = Y2 + shock2(i)
compute xi_2(i) = (Grossrate*Bond + G)/Y2_shocked(i)
compute C2(i)   = (1.0 - xi_2(i)-(1.0/(2.0*lambda))*(1 - sqrt(1 - 2.0*lambda*xi_2(i)))^2)*(Y2_shocked(i))
compute dU2(i)  = C2(i)^(-RiskAversion)*(1.0/sqrt(1.0-2.0*lambda*xi_2(i)))
compute E_dU2   = (1.0/(s_limit - lowerb))*%isimpson(Y2_shocked,dU2)
end do i

*******************
*** For figures ***
*******************
compute shock1_for_figure(j) = -shock1
compute Bond_for_figure(j) = %beta(1)

end find
end do j

*******************************************************************************
***                           Depicting a figure                            ***
*******************************************************************************
set shock_Fig 1 estimation_upper = shock1_for_figure(t - 2014:1 + 1)
set Bond_Fig 1 estimation_upper = Bond_for_figure(t - 2014:1 + 1)

scatter(style=line,key=below) 1
# shock_Fig Bond_Fig / 1