function [Lxx,Lxz,Llx,Llz,Lux,Luz,Clevel]=SolveLA(Cu,Cxl,Cz,Dxl,Fxl,Du,Fu,Dz,Fz,Pi,nx,nz) beta1=.95; alpha1=.33; W=-inv(Dxl-Du*inv(Cu)*Cxl)*(Fxl-Fu*inv(Cu)*Cxl); R=inv(Dxl-Du*inv(Cu)*Cxl)*((Dz+Du*inv(Cu)*Cz)*Pi+(Fz+Fu*inv(Cu)*Cz)); [T,S] = schur(W); eigW=eig(W); Tinv=inv(T); Q=Tinv*R; %u1=eigW(1); %u2=eigW(2); % q2=Q(2); % t21=Tinv(2,1); % t22=Tinv(2,2); % a1=(1/.95-((1-.33*.95)/(.33*.95))*(t21/t22)); % c1=1/(.33*.95); % c2=(1-.95*.33)/(.33*.95); % c3=q2/u2; % c4=t22*(1-Pi/u2); % a2=c1-c2*(c3/c4); % a2=((1/(.33*.95))-((1-.95*.33)/(.33*.95))*((q2/u2)/(t22*(1-Pi/u2)))); Txx=T(1:nx,1:nx); Sxx=S(1:nx,1:nx); Wxx=W(1:nx,1:nx); Wxl=W(1:nx,nx+1:size(W,2)); Rx=R(1:nx); TLL=Tinv(nx+1,nx+1); TLX=Tinv(nx+1,1:nx); Llx=-inv(TLL)*TLX; Llz=inv(TLL)*Pi; Lux=inv(Cu)*Cxl*[eye(nx);Llx]; Luz=inv(Cu)*Cxl*[zeros(nx,nz);Llz]+inv(Cu)*Cz; Lxx=Txx*Sxx*inv(Txx); Lxz=(Wxl*inv(TLL)*Pi+Rx); %simulate for a for a 1% increase in z z0=0; T=1000; t=1; k0=0; u=randn(T,1)*.01; ktbar=(alpha1*beta1).^(1/(1-alpha1)); ctbar=ktbar^(alpha1)-ktbar; kt0=ktbar; theta0=1; z01=1; B1=ctbar-Lux*ctbar-Luz*ctbar; B2=(ctbar/ktbar)*Lux; B3=(ctbar*Luz); while t