function ind = fisher(p,q,b,option)
%FISHER  computes the chain-weighted price or quantity index for base year b
%        ind=fisher(p,q,b,option) has TxN inputs of prices (p) and 
%        quantities (q) and returns, indices:
%
%                :
%             Ind[b-2]  = Ind[b-1]/Fisher[b-1]
%             Ind[b-1]  = Ind[b]/Fisher[b]
%             Ind[b]    = sum p[b]*q[b]         
%             Ind[b+1]  = Ind[b]*Fisher[b+1]
%             Ind[b+2]  = Ind[b+1]*Fisher[b+2]
%                :
%
%        where Fisher[t]=sqrt(Laspeyres[t]*Paasche[t]) and if option is set
%        to 1, then the quantity index is computed:
%                    
%             Laspeyres[t]  = sum p[t]*q[t]/sum p[t] q[t-1] 
%             Paasche[t]    = sum p[t-1]*q[t]/sum p[t-1] q[t-1] 
%
%        and, if the option is set to 2, the price index is computed:
%
%             Laspeyres[t]  = sum p[t]*q[t]/sum p[t-1] q[t] 
%             Paasche[t]    = sum p[t]*q[t-1]/sum p[t-1] q[t-1]
%

%        Ellen McGrattan, 7-2-15
%        Revised, erm, 7-2-15

[T0,N0]   = size(p);
[T,N]     = size(q);
if ((T0~=T) | (N0~=N));
  error('Check dimensions of price and quantity matrices') 
end;
if option==1;
  y       = q;
  x       = p;
else
  y       = p;
  x       = q;
end;

ycur      = y(2:T,:);
xcur      = x(2:T,:);
ylag      = y(1:T-1,:);
xlag      = x(1:T-1,:);
indl      = zeros(T,1);
indp      = zeros(T,1);

indl(2:T) = sum(xcur'.*ycur')'./sum(xcur'.*ylag')';
indp(2:T) = sum(xlag'.*ycur')'./sum(xlag'.*ylag')';
indf(2:T) = sqrt(indl(2:T).*indp(2:T));

ind(b)    = sum(x(b,:).*y(b,:));
for t=b+1:T;
  ind(t)  = ind(t-1)*indf(t);
end;
for t=b-1:-1:1;
  ind(t)  = ind(t+1)/indf(t+1); 
end;
ind       = ind(:);
