Inadequate matrix size for sensitivity analysis

Hi all,

I'm attempting to run an elasticity analysis on a matrix that's only 2 x 2:

syms AS JS Fe Pf
Svr = [AS JS Fe Pf];
mx = [0,   AS*Fe*Pf;
      JS,  AS];

However, when I run the code with this matrix design, the results do not come out right. When I run the code with a matrix that's 4 x 4, however, the code does work. A colleague told me that there's some issues with the zeroes, but she didn't know how to fix it. Any advice?

Here is the portion of the code that's not working right:

realmx=subs(mx,Svr,meanvr);
[lambdas,lambda1,W,w,V,v]=eigenall(realmx);
meanlam1=lambda1; sensmx=v*w'/(v'*w);
elastmx=(sensmx.*realmx)/lambda1;
meansens = zeros(1,numvrs);

for xx = 1: numvrs
    diffofvr=subs(diff(mx,Svr(xx)),Svr,meanvr);
    meansens(xx) = sum(sum(sensmx.*diffofvr));
end;
meanelast=((meansens.*meanvr)/lambda1);

maxlams=zeros(1,numvrs);
for rate = 1:numvrs
    vrates=meanvr; vrates(rate)=vrhi(rate);
    realmx=subs(mx,Svr,vrates);
    [lambdas,lambda1,W,w,V,v]=eigenall(realmx);
    maxlams(rate)=lambda1;
end;

disp(Svr); disp(meanelast);

for vr=1:numvrs
    x=sort(allelasts(:,vr));
    CLup(vr)=x(1+round((reps-1)*0.975));
    CLlo(vr)=x(1+round((reps-1)*0.025));
end;
disp(CLup); disp(CLlo);

The answer I get from that is:

meanelast = 0.8280    0.1720    0.1720    0.1720

CLup = 0.8291    0.3082    0.3082    0.3082
CLlo =    0.6918    0.1709    0.1709    0.1709

Which is erroneous because all elasticities should sum to 1. However, if I use a 4 x 4 matrix, I get elasticities that do sum to 1. Does anyone know how I might be able to fix this code so that the elasticities will process correctly with a 2 x 2 matrix?

Many, many thanks for taking the time to look at my issue!

[Star Strider]

What is ‘eigenall’?

I can’t find it in the MATLAB documentation. Could it be the problem?

[Sabrina]

Hi StarStrider,

Eigenall is the following code, which I run as a separate file:

[W,lambdas]=eig(mx);      

V=conj(inv(W));      

lambdas=diag(lambdas);    
[lambdas,I]=sort(lambdas);  

lambdas=flipud(lambdas);  

lambda1=lambdas(1);  
I=flipud(I);                
W=W(:,I);          
V=V(I,:);          
w=W(:,1);          
w=w/sum(w);          
v=real(V(1,:))';      
v=v/v(1);

Many thanks,

Sabrina