vnbrody.m

function [k, X, Y] = vnbrody(F, signoX, signoXF)
%vnbrody
%   [k, X, Y] = vnbrody(F)
%
%   Formatea el problema y lanza un c�digo de Br�dy para resolver VN.
%
%   [k, X, Y] = vnbrody(F, signoX, signoXF)
%      F         recetas en tiempo discreto
%      signoX    signo de las intensidades; por defecto X >= 0
%      signoXF   signo de los balances materiales; por defecto X F(k) = 0
%      k, X, Y   factor, intensidades y valores soluci�n
%
%   Los c�digos para los signos son   0 =    1 >=    2 <=    3 >=<
%   El algoritmo puede consultarse en:
%   Andr�s Br�dy, "The implicit dynamics of the Neumann growth model",
%   Acta Oeconomica, Vol. 54 (1) pp. 63-72 (2004)
%   http://www.akademiai.com/content/x71485711558/

%necesita formaproblema.m, formaab.m, formasolucion.m

if nargin < 2, signoX = []; end
if nargin < 3, signoXF = []; end

[Fn, signoX, signoXF] = formaproblema(F, signoX, signoXF, 1); %convierte el problema a la forma can�nica
[A, B] = formaab(Fn, 1);                                      %convierte a la forma A, B
[k, p, x] = neumann(A', B');                                  %lanza el c�digo original de Brody
X = x(1:size(Fn, 1));                                         %extrae las intensidades de los procesos que se inician
Y = p(1:size(Fn, 2))';                                        %extrae los precios de las materias habituales
[X, Y] = formasolucion(X, Y, signoX, signoXF, 1);             %convierte la soluci�n desde la forma can�nica





function [rho,p,x] = neumann(in,out)
% the function [rho,p,x] = neumann(in,out) computes a solution
% rho,p,x (growth factor, prices and quantities),
% for rectangular matrices in and out(input and output).

% Preliminaries:
% initial vectors, scaling, flag for indicating selection,
% initial rate, and select = 1 for going on with the process

[m,n] = size(out);
p = ones(1,m)/m;
x = ones(1,n)/n;

scale = m;
b = in/scale;
a = out/scale;

flag = [0 0 0 0 0];

rho = (p*a*x')/(p*b*x');
select = 1;

while select
   c = (a-rho*b);
   for i = 1:100      
      dx = p*c; %surplus gain
      dp = x*c'; %excess product
      p = p-dp;
      x = x+dx;
      
      p = p.*(p>0); %selecting valuable goods
      p = p/sum(p); %normalizing
      x = x.*(x>0); %selecting applied processes
      x = x/sum(x); %normalizing
      
      R(i,:) = p;
      Y(i,:) = x;
   end
   p = median(R); %taking median price
   x = median(Y); %taking median quantity
   p = p.*(p>0);
   p = p/sum(p);
   x = x.*(x>0);
   x = x/sum(x);
   rho = (p*a*x')/(p*b*x');
   
   %adjusting flag
   flag(1:4) = flag(2:5);
   flag(5) = sum(p>0) == sum(x>0); %indicating squareness
   
   %testing if the matrix was square for the last 5 iterations
   if sum(flag) == 5
      select = 0;
   end
end

%Selection finished. Finding equilbrium

pindex = find(p>0); %indexes of selected goods
xindex = find(x>0); %indexes of selected processes

A = a(pindex,xindex); %square matrix of outputs
B = b(pindex,xindex); %square matrix of inputs

%Inverting an almost singular matrix to improve exactness
for k = 1:3
   s = inv(A-rho*B);
   P = sum(s)/sum(s(:));
   X = sum(s')/sum(s(:));
   rho = P*A*X'/(P*B*X');
end

%Final values and exactness

p = zeros(1,m);
p(pindex) = P;
x = zeros(1,n);
x(xindex) = X;

sistem = rho*in-out;
error = num2str(p*sistem*x');
disp(['Exactness',error])