Setup

\(\(\( (*\ Copyright\ by\ Marco\ Bassetto, \ Roozbeh\ Hosseini, \ 
      Vadym\ Lepetyuk, \ 
      Christopher\ Phelan\ *) \)\(
\)\( (*\ \(This\ code\ can\ be\ freely\ distributed\ and\ modified\ for\ research\ purposes\ only\)\(,\)\ *) \)\(
\)\( (*\ 
      provided\ this\ copyright\ notice\ is\ included\ in\ the\ modified\ \(\(code\)\(.\)\)*) \)\(\ \)\(
\)\( (*\ 
      Proper\ credit\ should\ be\ given\ in\ all\ publications\ arising\ from\ *) \)\(
\)\( (*\ modifications\ of\ this\ code; \ 
      this\ should\ include\ a\ citation\ of\ *) \)\(
\)\( (*\ "Tax Riots"\ by\ Marco\ Bassetto\ and\ Christopher\ Phelan\ *) \)\)\)

Mathematica Setup

<< Graphics`Colors`
<< Graphics`Graphics`

Off[General::"spell"];
Off[General::"spell1"];
Off[Display::"pserr"];

If[$VersionNumber < 
        5.0, 
Total[x_]  := 
        Sum[x[[i]], {i, 1, Dimensions[x][[1]]}]]

SetOptions[Display, 
        ImageSize -> 360];

    SetOptions[Export, ImageSize -> 360];

$FormatType = TraditionalForm;

\(\(\( (*\ 
      current\ directory\ is\ set\ to\ match\ the\ location\ of\ this\ notebook; \ 
      eps\ files\ are\ saved\ there\ *) \)\(
\)\(\(filename = ToFileName["FileName" /. 
            NotebookInformation[EvaluationNotebook[]]];\)

    \(SetDirectory[DirectoryName[filename]];\)\)\)\)

Model Setup

u[c_ ]  :=  2 (Sqrt[c])  ;

\(\(\(\(pibar\ \  = \  .02;\)\ \ \  (*\ 
      maximum\ fraction\ of\ audited\ households\ *) 

    \(mlow\  = 0;\)\ \ \ \  (*\ 
      low\ bound\ of\ probability\ of\ low\ endowment\ *) 

    \(mhigh\  = 1;\)\)\(\ \ \ \ \)\( (*\ 
      upper\ bound\ of\ probability\ of\ low\ endowment\ *) \)\)\)

pi[m_]  :=  
        If[m > 0, Min[1. , (pibar/m)], 1];      (* 
      audit probability for a household which reports low endowment *)

elow  =  1/2 ;     
ehigh  =  
        3/2 ;

fpdf[m_] := 1 + 0  m

First Order Conditions

(* 
      first order condition with respect to tau1low *) 
eq1[tauhigh_, tau1low_, tau0low_, mu_, m_, 
        lambda_]  =  u'[elow  -  tau1low]  -  mu

(* 
      first order condition with respect to tauhigh *) 
eq3[tauhigh_, tau1low_, tau0low_, mu_, m_, 
        lambda_]  =  u'[ehigh - tauhigh] (1 + lambda)  -  
        mu

eq4[tauhigh_, tau1low_, 
        tau0low_, mu_, m_, lambda_]  =  
      m u'[elow - tau0low]  -  
        lambda (1 - m) u'[ehigh - tau0low]  -  
        m mu

eq4low[tauhigh_, tau1low_, 
        tau0low_, mu_, m_, lambda_]  =  
      mlow u'[elow - tau0low]  -  
        lambda (1 - mlow) u'[ehigh - tau0low]  -  
        mlow mu

eq4high[tauhigh_, tau1low_, 
        tau0low_, mu_, m_, lambda_]  =  
      mhigh  u'[elow - tau0low]  -  
        lambda (1 - mhigh) u'[ehigh - tau0low]  -  
        mhigh mu

eq5[tauhigh_, 
        tau1low_, tau0low_, mu_, m_, 
        lambda_]  =  (1 - m) tauhigh  +  
        m (pi[m] tau1low  +  (1 - pi[m]) tau0low)

eq5lowlow[tauhigh_, tau1low_, 
        tau0low_, mu_, m_, 
        lambda_]  =  (1 - 
              m) tauhigh  +  m tau1low

Private Information: Optimal Policy

Optimal policy over [plow;phigh]

Delta[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
          mustart_] := 
If[
          m > pibar, 
Deltamain[m, lambda, tauhighstart, 
            tau1lowstart, tau0lowstart, 
            mustart], 
Deltalowlow[m, lambda, tauhighstart, 
            tau1lowstart, mustart]
];

Optimal policy over [max{plow,pibar};phigh]

eqvec[tauhigh_, tau1low_, tau0low_, mu_, m_, 
          lambda_] = 
{eq1[tauhigh, tau1low, tau0low, mu, 
            m, lambda], 
eq3[tauhigh, tau1low, tau0low, mu, 
            m, lambda], 
eq4[tauhigh, tau1low, tau0low, mu, 
            m, lambda], 
eq5[tauhigh, tau1low, tau0low, mu, 
            m, lambda]};

policyrule[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
          mustart_] := 
FindRoot[
          eqvec[tauhigh, tau1low, tau0low, mu, m, 
            lambda], 
{tauhigh, tauhighstart}, {tau1low, 
            tau1lowstart}, {tau0low, tau0lowstart}, {mu, 
            mustart}];

Deltamain[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
        mustart_] := (pt = 
          policyrule[m, lambda, tauhighstart, tau1lowstart, tau0lowstart, 
            mustart]; 
u[ehigh - tauhigh] - (1 - pi[m]) 
              u[ehigh - tau0low] /. pt)

Optimal policy over [0;pibar]

eqveclowlow[tauhigh_, tau1low_, tau0low_, mu_, m_, 
          lambda_]  = 
{eq1[tauhigh, tau1low, tau0low, mu, 
            m, lambda], 
eq3[tauhigh, tau1low, tau0low, mu, 
            m, lambda], 
eq5lowlow[tauhigh, tau1low, 
            tau0low, mu, m, lambda]};

policyrulelowlow[m_, lambda_, tauhighstart_, tau1lowstart_, 
          mustart_]  := 
FindRoot[
          eqveclowlow[tauhigh, tau1low, 0, mu, m, 
            lambda], 
{tauhigh, tauhighstart}, {tau1low, 
            tau1lowstart}, {mu, mustart}];

Deltalowlow[m_, lambda_, tauhighstart_, tau1lowstart_, 
        mustart_]  :=  (pt  =  
          policyrulelowlow[m, lambda, tauhighstart, tau1lowstart, 
            mustart]; 
pt = 
          Join[pt, {tau0low -> 0}]; 
u[
            ehigh - tauhigh] /. pt)

Optimal policy over [pibar;plow]

eqveclow[tauhigh_, tau1low_, tau0low_, mu_, m_, 
          lambda_]  =  
{eq1[tauhigh, tau1low, tau0low, 
            mu, m, lambda], 
eq3[tauhigh, tau1low, tau0low, 
            mu, m, lambda], 
eq4low[tauhigh, tau1low, 
            tau0low, mu, m, lambda], 
eq5[tauhigh, tau1low, 
            tau0low, mu, m, lambda]};

policyrulelow[m_, lambda_, tauhighstart_, tau1lowstart_, 
          tau0lowstart_, mustart_]  :=  
        FindRoot[
          eqveclow[tauhigh, tau1low, tau0low, mu, m, 
            lambda], 
{tauhigh, tauhighstart}, {tau1low, 
            tau1lowstart}, {tau0low, tau0lowstart},  {mu, 
            mustart}];

Deltalow[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
        mustart_]  :=  (pt  =  
          policyrulelow[m, lambda, tauhighstart, tau1lowstart, tau0lowstart, 
            mustart]; 
u[
              ehigh - tauhigh]  -  (1 - pi[m]) u[ehigh - tau0low] /. 
          pt)

Optimal policy over [phigh;1]

eqvechigh[tauhigh_, tau1low_, tau0low_, mu_, m_, 
          lambda_]  = 
 {eq1[tauhigh, tau1low, tau0low, 
            mu, m, lambda], 
eq3[tauhigh, tau1low, tau0low, 
            mu, m, lambda], 
eq4high[tauhigh, tau1low, 
            tau0low, mu, m, lambda], 
eq5[tauhigh, tau1low, 
            tau0low, mu, m, lambda]};

policyrulehigh[m_, lambda_, tauhighstart_, tau1lowstart_, 
          tau0lowstart_, mustart_]  :=  
        FindRoot[
          eqvechigh[tauhigh, tau1low, tau0low, mu, m, 
            lambda], 
{tauhigh, tauhighstart}, {tau1low, 
            tau1lowstart}, {tau0low, tau0lowstart},  {mu, 
            mustart}];

Deltahigh[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
        mustart_]  :=  (pt  =  
          policyrulehigh[m, lambda, tauhighstart, tau1lowstart, tau0lowstart, 
            mustart]; 
u[
              ehigh - tauhigh]  -  (1 - pi[m]) u[ehigh - tau0low] /. 
          pt)

Private Information: Solution

Grid Setup

gm =  .0005;

Solution over grid on [plow;phigh]

mgrid = Table[m, {m, mlow, mhigh, gm}];
nm = Length[mgrid]

fvec = fpdf[mgrid];      
fvec = fvec/Total[fvec];

deltagrid = Table[0, {nm}];
tauhighgrid = Table[0, {nm}];
taulow1grid = Table[0, {nm}];
taulow0grid = Table[0, {nm}];
mugrid = Table[0, {nm}];

\(dist = 1. ; \ 
    adjcoef =  .5; \  (*\ 
      adjcoef\ might\ have\ to\ be\ adjusted\ depending\ on\ the\ parameters\ to\ achieve\ convergence\ *) 
tauhight0 =  .924; tau1lowt0 =  .1; 
    tau0lowt0 =  .1; 
    mut0 = 1.58;
\)
lambdat = 0;  

While[
      Abs[dist] >  .000001, 
(For[im = nm, 
          im ≥ 
            1, im--, 
(deltagrid[[im]] = 
              Delta[mgrid[[im]], lambdat, tauhight0, tau1lowt0, 
                tau0lowt0, 
                mut0]; 
tauhight0 = tauhighgrid[[im]] = tauhigh /. pt; 
tau1lowt0 = taulow1grid[[im]] = 
                tau1low /. 
                  pt; 
tau0lowt0 = taulow0grid[[im]] = 
                tau0low /. 
                  pt; 
mut0 = mugrid[[im]] = 
                mu /. pt)]; 
tauhight0 = 
          tauhighgrid[[nm]];  tau1lowt0 = taulow1grid[[nm]];  
        tau0lowt0 = taulow0grid[[nm]];  
        mut0 = mugrid[[nm]]; 
weighteddelta = 
          deltagrid (1 - mgrid) fvec; 
dist = 
          Total[weighteddelta]; 
lambdat = 
          lambdat - adjcoef dist; 
Print[{dist, 
            lambdat}];)]

Solution over grid on (pibar;plow)

mgridlow  =  
        Table[m, {m, pibar + gm, mlow - gm, 
            gm}];
nmlow  =  Length[mgridlow]

deltalowgrid  =  Table[0, {nmlow}];
tauhighlowgrid  =  Table[0, {nmlow}];
taulow1lowgrid  =  Table[0, {nmlow}];
taulow0lowgrid  =   Table[0, {nmlow}];
mulowgrid  =  Table[0, {m, nmlow}];

tauhight0 =  tauhighgrid[[1]];
tau1lowt0 =   taulow1grid[[1]];
tau0lowt0  =  taulow0grid[[1]];
mut0 =  mugrid[[1]];

For[im = nmlow, 
        im ≥ 
          1, im--, 
(deltalowgrid[[im]]  =  
            Deltalow[mgridlow[[im]], lambdat, tauhight0, tau1lowt0, 
              tau0lowt0,  
              mut0]; 
tauhight0 = tauhighlowgrid[[im]] = 
              tauhigh /. 
                pt; 
tau1lowt0 = taulow1lowgrid[[im]] = 
              tau1low /. 
                pt; 
tau0lowt0 = taulow0lowgrid[[im]] = 
              tau0low /. 
                pt; 
mut0 = mulowgrid[[im]] = 
              mu /. pt)];

Solution over grid on [0;min{pibar,plow}]

mgridlowlow  =  
        Table[m, {m, 0, Min[pibar, mlow - gm], 
            gm}];
nmlowlow  =  Length[mgridlowlow]

deltalowlowgrid  =  Table[0, {nmlowlow}];
tauhighlowlowgrid  =  Table[0, {nmlowlow}];
taulow1lowlowgrid  =  Table[0, {nmlowlow}];
taulow0lowlowgrid  =   Table[0, {nmlowlow}];
mulowlowgrid  =  Table[0, {nmlowlow}];

tauhight0 =  tauhighgrid[[1]];
tau1lowt0 =   taulow1grid[[1]];
mut0 =  mugrid[[1]];

For[im = nmlowlow, 
        im ≥ 
          1, im--, 
(deltalowlowgrid[[im]]  =  
            Deltalowlow[mgridlowlow[[im]], lambdat, tauhight0, tau1lowt0, 
              mut0]; 
tauhight0  =  tauhighlowlowgrid[[im]] = 
              tauhigh /. 
                pt; 
tau1lowt0  =  taulow1lowlowgrid[[im]] = tau1low /. 
                pt; 
mut0  =  mulowlowgrid[[im]] = mu /. pt)];

Solution over grid on (phigh;1]

mgridhigh  =  
        Table[m, {m, mhigh + gm, 1, gm}];
nmhigh  =  Length[mgridhigh]

deltahighgrid  =  Table[0, {nmhigh}];
tauhighhighgrid  =  Table[0, {nmhigh}];
taulow1highgrid  =  Table[0, {nmhigh}];
taulow0highgrid  =   Table[0, {nmhigh}];
muhighgrid  =  Table[0, {nmhigh}];

tauhight0 =  tauhighgrid[[nm]];
tau1lowt0 =   taulow1grid[[nm]];
tau0lowt0  =  taulow0grid[[nm]];
mut0 =  mugrid[[nm]];

For[im  =  1, 
        im ≤ 
          nmhigh, im++, 
(deltahighgrid[[im]]  =  Deltahigh[mgridhigh[[im]], lambdat, tauhight0, tau1lowt0, tau0lowt0, 
              mut0]; 
tauhight0 = tauhighhighgrid[[im]] = 
              tauhigh /. 
                pt; 
tau1lowt0 = taulow1highgrid[[im]] = tau1low /. 
                pt; 
tau0lowt0 = taulow0highgrid[[im]] = tau0low /. pt; 
mut0 = muhighgrid[[im]] = 
              mu /. pt)];

Join Grid

nmfull  =  nmlowlow + nmlow + nm + nmhigh

mgridfull  =  
        Join[mgridlowlow, mgridlow, mgrid, 
          mgridhigh];
deltagridfull   =  
        Join[deltalowlowgrid, deltalowgrid, deltagrid, 
          deltahighgrid];
tauhighgridfull   =  
        Join[tauhighlowlowgrid, tauhighlowgrid, tauhighgrid, 
          tauhighhighgrid];
taulow1gridfull   =  
        Join[taulow1lowlowgrid, taulow1lowgrid, taulow1grid, 
          taulow1highgrid];
taulow0gridfull   =  
        Join[taulow0lowlowgrid, taulow0lowgrid, taulow0grid, 
          taulow0highgrid];
mugridfull   =  
        Join[mulowlowgrid, mulowgrid, mugrid, 
          muhighgrid];
fvecfull  =  
        Join[Table[0, {nmlowlow + nmlow}], fvec, 
          Table[0, {nmhigh}]];

Private Information: Graphs

pstyle = {Red, AbsoluteThickness[1.6]} ;

deltag  =  
        ListPlot[
          Transpose[{mgridfull, 
              deltagridfull}], 
PlotJoined -> True, 
          PlotRange -> All, PlotStyle -> pstyle, 
          Epilog -> {Text["m", { .9,  .15}], 
              Text["△(m)", {- .15, 
                  1}]}];
Export["deltag.eps", deltag, "eps"];

fg  =  
        ListPlot[
          Table[{mgridfull[[im]], fvecfull[[im]]}, {im, 1, 
              nmfull}], 
PlotJoined -> True, 
          PlotRange -> All, 
          PlotStyle -> {Cobalt, Thickness[ .006]}];

tauhighg  =  
        ListPlot[
          Transpose[{mgridfull, 
              tauhighgridfull}], 
PlotJoined -> True, 
          AxesOrigin -> {0, 0}, 
          PlotRange -> {{0, 1}, {- .05, 1}}, 
          PlotStyle -> pstyle];

taulow1g  =  
        ListPlot[
          Transpose[{mgridfull, 
              taulow1gridfull}], 
PlotJoined -> True, 
          PlotRange -> All, PlotStyle -> pstyle];

ipibar = Ceiling[pibar/gm] + 1;
taulow0g  =  
        ListPlot[
          Table[{mgridfull[[im]], taulow0gridfull[[im]]}, {im, 
              ipibar + 1, nmfull}], 
PlotJoined -> 
            True, PlotRange -> All, AxesOrigin -> {0, 0}, 
          PlotStyle -> pstyle];

Private Information: Effects of redistribution

Grid Setup

imlow = Ceiling[mlow/gm] + 1;    
imlow = Max[imlow, 2];      
nmrd = 
      nmfull - imlow;

xgrid = Table[0, {nmrd}];
incentivegrid = Table[0, {nmrd}];

Incentives under redistribution

For[im = imlow, 
        im ≤ 
          nmfull - 
            1, im++, 
(mtildalow = (im - 1) 
              gm; 
ix = 
            im - imlow + 
              1; 
xgrid[[ix]] = x = (mtildalow - 
                    mlow)/(1 - mlow); 
mtildahigh = 
            mhigh + (1 - mhigh) x; 
imtildahigh = 
            Min[Floor[mtildahigh/gm] + 1, 
              nmfull]; 
incentivegrid[[ix]] = 
            Sum[
(u[
                        ehigh - tauhighgridfull[[iim]] + 
                          pibar x/(1 - x)/
                              mgridfull[[iim]] (ehigh - 
                                taulow1gridfull[[iim]])] - (1 - 
                            pi[mgridfull[[iim]]]) 
                        u[ehigh - taulow0gridfull[[iim]]]) (1 - 
                      mgridfull[[iim]]) 
                  fpdf[(mgridfull[[iim]] - x)/(1 - x)], {iim, im, 
                  imtildahigh}]/(1 - 
                    x)^2;
)];

ListPlot[
        Transpose[{xgrid, incentivegrid}], {PlotJoined -> True, 
          PlotStyle -> Blue}];

(* 
      Searching for x such that incentive[x] = 
        0 *) 
For[ix = nmrd, 
      ix ≥ 1, ix--, 
      If[incentivegrid[[ix]] < 
          0, 
(Print[{xgrid[[ix]], 
              incentivegrid[[ix]]}]; 
Print[{xgrid[[ix + 1]], 
              incentivegrid[[ix + 
                    1]]}]; 
Break[])]]

\(\(\( (*\ 
      linear\ approximation\ *) \)\(
\)\(\((incentivegrid[\([ix]\)] 
            xgrid[\([ix + 1]\)] - 
          incentivegrid[\([ix + 1]\)] 
            xgrid[\([ix]\)])\)/\((incentivegrid[\([ix]\)] - 
          incentivegrid[\([ix + 1]\)])\)\)\)\)