clear all
clc

N = 2 ; Y = [ 1  4 ]'; 
e = 5 ;

lambda = [0.88 0.02 ; ... 
          0.12 0.98];

tau_hat = [1 0.8 ; ...
           0.8 1];

 x0 = ones(N,1);
 w_hat = fsolve(@(x) HAT_ALGEBRA(x, tau_hat, lambda , N, Y, e), x0);
 
 
 c_hat = tau_hat.* repmat(w_hat,1,N);
 lambda_new = lambda.*( c_hat.^-e )./repmat(sum(lambda.* (c_hat.^-e) ),N,1);
 W_hat = ( lambda_new(eye(N)==1)./lambda(eye(N)==1) ).^(-1/e);
 Welfare_Gain_pct = 100 .* (  W_hat - 1 );
 disp(['change in US Real GDP p/c = ', num2str(Welfare_Gain_pct(1)), '%'])
 disp(['change in ROW Real GDP p/c = ', num2str(Welfare_Gain_pct(2)), '%'])

  function ERR = HAT_ALGEBRA(x, tau_hat, lambda , N, Y, e)
    c_hat = tau_hat .* repmat(x,1,N);
    SALES = lambda .* (c_hat.^-e) ./ repmat( sum(lambda.* (c_hat.^-e) ), N,1) .* repmat(x' .* Y',N,1);
    ERR =  x.* Y - sum( SALES , 2);
  end     
  
  
  
