function burn
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% computes potential money burning allocation for 3 type case
% reports allocation and whether or not it satisfies
% auxiliary conditions to be deemed IC
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%% parameterization 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sigma=2, thetal = .8, thetah=1.6, thetam=1, ploverph = .7 betta=.69
% leads to money burning for the middle type for high enough pm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
global betta sigma thetal thetah thetam lamda y p

% parameters
ploverph = .7;  %ratio of pl over ph

ppmm=(.001:.001:.999)';  %possible values for pm
sigma = 2;
thetal = .8;
thetah = 1.6;
thetam = 1;
betta = .7

y=1;


% first best allocation
cl_fb = y*((thetal)^(-1/sigma)+1)^(-1);
ch_fb = y*((thetah)^(-1/sigma)+1)^(-1);
kl_fb = y - cl_fb;
kh_fb = y - ch_fb;

beta_star = thetal*(u(cl_fb) - u(ch_fb) ) / (u(kh_fb) - u(kl_fb));
% when betta<beta_star, first best for low and high not IC
Ym =0;  % dummy variable initialized


if (betta > thetal/thetam | betta < thetam/thetah )
    'no money burning possible-!!!! change parameters'
else

    for i = 1: length(ppmm)
        %    pause
        clc;
        pm=ppmm(i);

        %construction ph and pl given pm and the ratio pl/pm
        ph=(1-pm)/(1+ploverph);
        pl=1-pm-ph;


        % solve for mu (eq 3-4)
        A = [ -thetal/thetam, 1 ;  - betta , betta];
        B = [ -pm               ; -pm];
        mu=inv(A)*B;
        if mu(1) < 0 | mu(2) < 0; display('mu1 or mu2 is negative'); problem(i,1)=1; else

            % compute cL and cH from mu (eq 5-6)
            rl = thetal*(pl + mu(1))/(pl + betta*mu(1));  % ratio of marginal utilities for the low type
            rh = thetah*(ph - mu(2)*thetam/thetah)/(ph - betta*mu(2) ); % ratio of marginal utilities for the high type
            % if any of these ratios is negative .. this cannot be possible.. stop
            if rl<0 | rh<0; display('rl or rh are negative'); problem(i,1)=1.5; else

                %compute the allocation for low and high from the ratio of marginal utilities and income
                c1l = y*[ 1 + rl.^(-1/sigma)].^-1; c2l = y - c1l;
                c1h = y*[ 1 + rh.^(-1/sigma)].^-1; c2h = y - c1h;

                u1l = u(c1l); u2l = u(c2l); u1h = u(c1h); u2h = u(c2h);

                % now compute u1m and u2m from binding linear IC equations for low and medium types (
                % (equations 7-8)
                Aic = [ thetal , betta ;  thetam , betta];
                Bic = [ thetal*u1l + betta*u2l ; thetam*u1h + betta*u2h];

                um = inv(Aic)*Bic;

                %checking that those utility values are feasible
                if (sigma-1)*um(1)>0 | (sigma-1)*um(2) >0
                    ' utility is out of bounds'
                    problem(i,1) = 2;
                else

                    %finding the consumption bundle for the middle type
                    c1m =[(1-sigma)*um(1)].^(1/(1-sigma));
                    c2m =[(1-sigma)*um(2)].^(1/(1-sigma));
                    c1  = [c1l, c1m, c1h]
                    C1(i,:) = c1;

                    % check monotonicity, which is a necessary condition for IC and has not been imposed yet
                    if c1l > c1h | c1m > c1h | c1l > c1m | c2m > c2l | c2m < c2h | c2l < c2h
                        'proposed solution is not IC!'
                        problem(i,1)=3;
                    else

                        % check sign of lamdaH (multiplier of the resource for the high type)
                        %             the multiplier for the low type is positive if mu is positive
                        if ph - thetam/thetah*mu(2) < 0 | ph - betta*mu(2) < 0 ;
                            'lamdam turned out negative'
                            problem(i,1) = 4;
                        else


                            ym = c1m + c2m

                            Ym(i,1) = ym;
                            u1m = u(c1m); u2m = u(c2m);
                            U = pl*[thetal*u1l + u2l] + pm*[thetam*u1m + u2m]+ ph*[thetah*u1h + u2h];
                            Etheta = (pl*thetal+pm*thetam+ph*thetah);
                            cpool = (1+Etheta^(-1/sigma))^(-1);
                            Upool = Etheta*u(cpool)+u(1-cpool);

                            if Upool > U
                                'pooling is better :p'
                                problem(i,1) = 5;
                            else
                                'separating is better than pooling'

                            end
                        end
                    end
                end
            end
        end
    end
end

clf
%subplot(2,1,1)
figure(1)
titletext(1) = {'c1 allocation for l, m and h '};
titletext(2) = {['\sigma =',num2str(sigma),' \theta_l =',num2str(thetal),' \thetah =',num2str(thetam),' pl/ph =',num2str(ploverph),' betta =',num2str(betta)]};
beta_star

plot(ppmm(find(problem==0)) , C1(find(problem==0),:))
grid
xlabel('pm')
ylabel('c_1')
legend('c_1_l','c_1_m','c_1_h')
title( titletext , 'fontweight' , 'bold')

figure(2)
%subplot(2,1,2)
plot(ppmm(find(problem==0)) , Ym(find(problem==0)))
grid
xlabel('pm')
ylabel('c_1_m + c_2_m ')
title('total consumption for m-type','fontweight','bold')




function f=u(x)

global sigma

f = (1/(1-sigma))*x.^(1-sigma);