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\ *) \)\)\)

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

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

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

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

    \(mbar\  =  .5;\)\)\(\ \ \ \ \)\( (*\ 
      probability\ of\ low\ endowment\ *) \)\)\)

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

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

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

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

First Order Conditions

eq1[tauhigh_, tau1low_, tau0low_, mu_, m_, 
        lambda_]  =  u'[elow  -  tau1low]  -  mu

eq3[tauhigh_, tau1low_, tau0low_, mu_, m_, 
        lambda_]  =  u'[ehigh - tauhigh] (1 + lambda)  -  
        mu

eq4[tauhigh_, tau1low_, tau0low_, mu_, m_, lambda_]  =  
      mbar u'[elow - tau0low]  -  
        lambda (1 - mbar) u'[ehigh - tau0low]  -  
        mbar 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

Optimal Policy

Optimal policy over [pibar;1]

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]};

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

policyvec[1,  .1,  .9,  .1,  .01, 1.5]

Delta[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
        mustart_] := (pt = 
          policyvec[m, lambda, tauhighstart, tau1lowstart, tau0lowstart, 
            mustart]; 
tauhight = pt[[1]]; 
        tau0lowt = 
          pt[[3]]; 
u[
            ehigh - tauhight] - (1 - pi[m]) 
            u[ehigh - tau0lowt])

Delta[1,  .2,  .924,  .1,  .01, 1.58114]

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]};

policyveclowlow[m_, lambda_, tauhighstart_, tau1lowstart_, 
        tau0lowstart_, 
        mustart_]  :=  (policyrule  =  
          FindRoot[
            eqveclowlow[tauhigh, tau1low, 0. , mu, m, 
              lambda], 
{tauhigh, tauhighstart}, {tau1low, 
              tau1lowstart}, {mu, mustart}]; 
{tauhigh, 
            tau1low, 0. , mu} /. policyrule)

Deltalowlow[m_, lambda_, tauhighstart_, tau1lowstart_, tau0lowstart_, 
        mustart_]  :=  (pt  =  
          policyveclowlow[m, lambda, tauhighstart, tau1lowstart, 
            tau0lowstart, mustart]; 
tauhight  = 
          pt[[1]]; 
u[ehigh - tauhight] )

Private Information: Solution

Solution at mbar

dist = 1. ;
tauhight0 =  .924; tau1lowt0 =  .1; tau0lowt0 =  .1; 
    mut0 = 1.58;
lambdat =  .30947;
While[
      Abs[dist] >  .000001, 
(deltat = 
          Delta[mbar, lambdat, tauhight0, tau1lowt0, tau0lowt0, 
            mut0]; 
tauhight0 = pt[[1]]; 
        tau1lowt0 = pt[[2]]; tau0lowt0 = pt[[3]]; 
        mut0 = pt[[4]]; 
dist = 
          deltat; 
Print[
          dist]; 
lambdat = lambdat - dist)]
tauhighmbar = pt[[1]]; tau1lowmbar = pt[[2]]; 
    tau0lowmbar = pt[[3]]; mumbar = pt[[4]];

lambdat

Solution over grid on [pibar;1]

gm =  .0005;

mgrid  =  Table[m, {m, pibar + gm, 1, gm}]; 
    nm  =  Length[mgrid]

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

tauhight0 =  tauhighmbar;
tau1lowt0 =   tau1lowmbar;
tau0lowt0  =  tau0lowmbar;
mut0 =  mumbar;

For[im  =  nm, 
        im ≥ 
          1, im--, 
(deltagrid[[im]]  =  
            Delta[mgrid[[im]], lambdat, tauhight0, tau1lowt0, tau0lowt0, 
              mut0]; 
tauhight0  =  pt[[1]]; 
          tau1lowt0  =  pt[[2]];   tau0lowt0  =  pt[[3]];  
          mut0  =  
            pt[[4]]; 
tauhighgrid[[im]] =  
            pt[[1]];  taulow1grid[[im]]  =  pt[[2]]; 
          taulow0grid[[im]]  =  
            pt[[3]]; 
mugrid[[im]]  =  
            pt[[4]])];

Solution over grid on [0;pibar]

mgridlowlow  =  Table[m, {m, 0, pibar, 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, 
              0. , mut0]; 
tauhight0  =  pt[[1]]; 
          tau1lowt0  =  pt[[2]];   mut0  =  
            pt[[4]]; 
tauhighlowlowgrid[[im]] =  
            pt[[1]];  
          taulow1lowlowgrid[[im]]  =  
            pt[[2]]; 
mulowlowgrid[[im]]  =  
            pt[[4]])];

Join Grid

nmfull  =  nmlowlow + nm

mgridfull  =  
        Join[mgridlowlow, mgrid];
deltagridfull   =  
        Join[deltalowlowgrid, deltagrid];
tauhighgridfull   =  
        Join[tauhighlowlowgrid, tauhighgrid];
taulow1gridfull   =  
        Join[taulow1lowlowgrid, taulow1grid];
taulow0gridfull   =  
        Join[taulow0lowlowgrid, taulow0grid];
mugridfull   =  Join[mulowlowgrid, mugrid];

Private Information: Graphs

deltag  =  
        ListPlot[
          Table[{mgridfull[[im]], deltagridfull[[im]]}, {im, 1, 
              nmfull}], 
PlotJoined -> True, 
          PlotRange -> All, 
          PlotStyle -> {Red, Thickness[ .01]}];

Display["deltag0.eps", 
        Show[deltag, 
          Graphics[Text["m", { .9, - .2}]]], "eps"];

taulow1g  =  
        ListPlot[
          Table[{mgridfull[[im]], taulow1gridfull[[im]]}, {im, 1, 
              nmfull}], 
PlotJoined -> True, 
          PlotRange -> {{0, 1}, {-0.4, 0}}, 
          PlotStyle -> {Red, Thickness[ .01]}, 
          AxesOrigin -> {0, 0}];

tauhighg  =  
        ListPlot[
          Table[{mgridfull[[im]], tauhighgridfull[[im]]}, {im, 1, 
              nmfull}], 
PlotRange -> All, 
          PlotJoined -> True, PlotRange -> All, 
          PlotStyle -> {Red, Thickness[ .01]}, 
          AxesOrigin -> {0, 0}];

Display["tauhighg0.eps", 
        Show[tauhighg, 
          Graphics[Text["m", { .5,  .1}]]], "eps"];

taulow0g  =  
        ListPlot[
          Table[{mgridfull[[im]], taulow0gridfull[[im]]}, {im, 
              nmlowlow + 1, nmfull}], 
PlotRange -> 
            All, PlotJoined -> True, PlotRange -> All, 
          PlotStyle -> {Red, Thickness[ .01]}];

Full Information: Solution

Solution over grid on [pibar;1]

deltagridf  =  Table[0, {nm}];
tauhighgridf  =  Table[0, {nm}];
taulow1gridf  =  Table[0, {nm}];
taulow0gridf  =  Table[0, {nm}];
mugridf  =  Table[0, {nm}];

lambdat  =  0. ;

tauhight0  =   .924;  tau1lowt0 =  .1; tau0lowt0  =   .1;  
    mut0  =  1.58;

For[im = nm, 
        im ≥ 
          1, im--, 
(deltagridf[[im]]  =  
            Delta[mgrid[[im]], lambdat, tauhight0, tau1lowt0, tau0lowt0, 
              mut0]; 
tauhight0  =  pt[[1]]; 
          tau1lowt0  =  pt[[2]];  tau0lowt0  =  pt[[3]];  
          mut0  =  
            pt[[4]]; 
tauhighgridf[[im]] =  
            pt[[1]];  
          taulow1gridf[[im]]  =  
            pt[[2]]; 
taulow0gridf[[im]] = 
            pt[[3]]; 
          mugridf[[im]]  =  pt[[4]])];

Solution over grid on [0;pibar]

deltalowlowgridf  =  Table[0, {nmlowlow}];
tauhighlowlowgridf  =  Table[0, {nmlowlow}];
taulow1lowlowgridf  =  Table[0, {nmlowlow}];
taulow0lowlowgridf  =   Table[0, {nmlowlow}];
mulowlowgridf  =  Table[0, {nmlowlow}];

tauhight0 =  tauhighgridf[[1]];
tau1lowt0 =   taulow1gridf[[1]];
mut0 =  mugridf[[1]];

For[im = nmlowlow, 
        im ≥ 
          1, im--, 
(deltalowlowgridf[[im]]  =  
            Deltalowlow[mgridlowlow[[im]], lambdat, tauhight0, tau1lowt0, 
              0. , mut0]; 
tauhight0  =  pt[[1]]; 
          tau1lowt0  =  pt[[2]];    mut0  =  
            pt[[4]]; 
tauhighlowlowgridf[[im]] =  
            pt[[1]];  
          taulow1lowlowgridf[[im]]  =  
            pt[[2]]; 
mulowlowgridf[[im]]  =  
            pt[[4]])];

Join Grid

deltagridfullf   =  
        Join[deltalowlowgridf, deltagridf];
tauhighgridfullf   =  
        Join[tauhighlowlowgridf, tauhighgridf];
taulow1gridfullf   =  
        Join[taulow1lowlowgridf, taulow1gridf];
taulow0gridfullf   =  
        Join[taulow0lowlowgridf, taulow0gridf];
mugridfullf   =  Join[mulowlowgridf, mugridf];

Full Information: Graphs

deltafg  =  
        ListPlot[
          Table[{mgridfull[[im]], deltagridfullf[[im]]}, {im, 1, 
              nmfull}], 
PlotJoined -> True, 
          PlotRange -> All, 
          PlotStyle -> {Blue, Thickness[ .01]}];

Show[deltag, deltafg];

taulow1fg  =  
        ListPlot[
          Table[{mgridfull[[im]], taulow1gridfullf[[im]]}, {im, 1, 
              nmfull}], 
PlotJoined -> True, 
          PlotRange -> All, PlotStyle -> {Blue, Thickness[ .01]}, 
          PlotRange -> {{0, 1}, {- .7,  .5}}];

Display["taulow1fg0.eps", 
        Show[taulow1fg, Graphics[Text["m", { .5,  .075}]], 
          PlotRange -> All], "eps"];

Display["taulow1g0.eps", 
        Show[taulow1fg, taulow1g, Graphics[Text["m", { .5,  .075}]], 
          PlotRange -> All], "eps"];

tauhighfg  =  
        ListPlot[
          Table[{mgridfull[[im]], tauhighgridfullf[[im]]}, {im, 1, 
              nmfull}], 
PlotRange -> All, 
          PlotJoined -> True, PlotRange -> All, 
          PlotStyle -> {Blue, Thickness[ .01]}];

Display["tauhighfg0.eps", 
        Show[tauhighfg, Graphics[Text["m", { .5,  .1}]], 
          AxesOrigin ->  {0, 0}], "eps"];

Display["tauhighg0.eps", 
        Show[tauhighfg, tauhighg, Graphics[Text["m", { .5,  .1}]], 
          AxesOrigin ->  {0, 0}], "eps"];

taulow0fg  =  
        ListPlot[
          Table[{mgridfull[[im]], taulow0gridfullf[[im]]}, {im, 
              nmlowlow + 1, nmfull}], 
PlotJoined -> 
            True, PlotRange -> {{0, 1}, {-0.3, 0.05}}, 
          PlotStyle -> {Blue, Thickness[ .01]}];

Display["taulow0fg0.eps", 
        Show[taulow0fg, 
          Graphics[Text["m", { .5, - .1}]]], "eps"];

Display["taulow0g0.eps", 
        Show[taulow0fg, taulow0g, 
          Graphics[Text["m", { .5, - .1}]]], "eps"];

Perturbing the solution on small-probability events

Cdfn[x_] = (1 + Erf[x/Sqrt[2]])/2;
InverseCDFn[p_] = Sqrt[2] InverseErf[0, 2  p - 1];

npeople = 1000000

We set  'mhat'  so that under truthtelling the probability of  'm>mhat'  is 1%.

mhat = 
      mbar + InverseCDFn[1 - 0.01] 
          Sqrt[mbar (1 - mbar)/npeople]

deltaalt = deltagridfull;
tauhighalt = tauhighgridfull;
taulow0alt = taulow0gridfull;

For[im = 1, 
        im ≤ nmfull - 1, im++, 
If[
          mgridfull[[im]] ≥ 
            mhat, 
deltaalt[[im]] = 
            u[ehigh + 
                mgridfull[[im]] elow/(1 - mgridfull[[im]])];  
          tauhighalt[[im]] = -mgridfull[[im]] 
              elow/(1 - mgridfull[[im]]);  
          taulow0alt[[im]] = elow]];

deltaalt[[nmfull]] = deltaalt[[nmfull - 1]]

We will use a linear approximation of the deltas around the relevant points. This is OK as long as x remains fairly small.

igrid = Floor[mhat/gm] + 1;
ialt = igrid + 1;

deltagridbase = deltaalt[[igrid]]
deltaaltbase = deltaalt[[ialt]]

deltagridinc = 
      deltaalt[[igrid]] - deltaalt[[igrid - 1]]
deltaaltinc = 
      deltaalt[[ialt + 1]] - deltaalt[[ialt]]

deltaavg = Table[0, {ix, 101}];

For[ix = 0, 
        ix ≤ 100, ix++, 
x = 
          ix/10000; μ = mbar + x (1 - mbar);  σ = 
          Sqrt[μ (1 - μ)/
                npeople]; 
deltaavg[[ix + 
                1]] = 
(deltagridbase + 
                  deltagridinc (μ - mgridfull[[igrid]])) 
              Cdfn[(mhat - μ)/σ]
 - (deltagridinc σ/Sqrt[2  π]) 
              Exp[-((mhat - μ)/σ)^2/
                  2]
 + (deltaaltbase + 
                  deltaaltinc (μ - mgridfull[[ialt]])) (1 - 
                  Cdfn[(mhat - μ)/σ])
 + (deltaaltinc σ/Sqrt[2  π]) 
              Exp[-((mhat - μ)/σ)^2/2]];

deltaavg

Perturbed solution: Graphs

deltax  =  
        ListPlot[
          Table[{ix/10000, deltaavg[[ix + 1]]}, {ix, 1, 
              100}], 
PlotJoined -> True, 
          PlotStyle -> {Red, Thickness[ .01]}, 
          PlotRange -> {{0,  .01}, {-1, 5}}];

Display["deltax0.eps", 
        Show[deltax, 
          Graphics[
            Text["x", { .009, - .2}]]], "eps"];