Merge branch 'NoImsl' of git@github.com:frbatlanta/TZcode into NoImsl

MatlabFiles/MSV/fwz_msv_msre.m

@@ -17,8 +17,10 @@ function [F1, F2, G1, G2, V, err] = fwz_msv_msre(P, A, B, Psi, s, x, max_count,
 %
 %  eta(t) = F2{s(t)}*x(t-1) + G2{s(t)}*epsilont(t)
 %
-% err is set to 0 upon sucess.  A positive value is the number of
-% iterations before the method was terminated.
+% A positive value of err is the number of iterations needed to obtain 
+% convergence and indicates success.  A negitive value of err is the number of 
+% iterations before the method terminated without convergence and indicates
+% failure.
 %
 
 h=size(P,1);
@@ -38,20 +40,9 @@ if nargin <= 5
 end
 f=zeros(h*s*(n-s),1);
 
-
 Is=eye(s);
 I=eye(n-s);
 
-% The following code would be used if Pi were passed instead of the
-% assumption that Pi' = [zeros(s,n-s)  eye(s)].
-% for i=1:h
-%     [Q,R] = qr(Pi{i});
-%     W=[zeros(n-s,s)  I; inv(R(1:s,1:s)); zeros(s,n-s)]*Q';
-%     A{i}=W*A{i};
-%     B{i}=W*B{i};
-%     Psi{i}=W*Psi{i};
-% end
-
 U=cell(h,1);
 for i=1:h
     U{i}=inv(A{i});
@@ -109,9 +100,9 @@ while cont
 end
 
 if (norm(f) < tol)
-    err=0;
-else
     err=count;
+else
+    err=-count;
 end
 
 F1=cell(h,1);

MatlabFiles/MSV/fwz_pi2eye.m

+function [A_new, B_new, Psi_new] = fwz_pi2eye(A, B, Psi, Pi)
+% [A, B, Psi] = fwz_pi2eye(A, B, Psi, Pi)
+% Converts the representation 
+%
+%   A(s(t) x(t) = B(s(t) x(t-1) + Psi(s(t)) epsilon(t) + Pi{s(t)} eta(t)
+%
+% to a form compatible with fwz_msv_msre().  This is the form
+%
+%   A(s(t) x(t) = B(s(t) x(t-1) + Psi(s(t)) epsilon(t) + Pi eta(t)
+%
+% where Pi' = [zeros(s,n-s)  eye(s)]
+%
+
+for i=1:h
+    [Q,R] = qr(Pi{i});
+    W=[zeros(n-s,s)  I; inv(R(1:s,1:s)); zeros(s,n-s)]*Q';
+    A_new{i}=W*A{i};
+    B_new{i}=W*B{i};
+    Psi_new{i}=W*Psi{i};
+end
\ No newline at end of file

MatlabFiles/MSV/verify_solution.m

+function diff = verify_solution(P,A,B,Psi,s,F1,F2,G1,G2,V)
+%
+% Verifies that 
+%
+%     x(t) = V(s(t))(F1(s(t))x(t-1) + G1(s(t))epsilon(t)
+%   eta(t) = F2(s(t))x(t-1) + G2(s(t))epsilon(t)
+%
+% is a solution of 
+% 
+%    A(s(t)) x(t) = B(s(t)) x(t-1) + Psi(s(t) epsilon(t) + Pi eta(t)
+%
+% where Pi' = [zeros(s,n-s)  eye(s)] by verifying that
+%
+%  [A(j)*V(j)  Pi] [ F1(j)   
+%                    F2(j) ] = B(j)
+%
+%  [A(j)*V(j)  Pi] [ G1(j)   
+%                    G2(j) ] = Psi(j)
+%
+%  (P(i,1)*F2(1) + ... + P(i,h)*F2(h))*V(i) = 0
+%
+
+h=size(P,1);
+n=size(A{1},1);
+r=size(Psi{1},2);
+
+diff=0;
+Pi=[zeros(n-s,s); eye(s)];
+for j=1:h
+    X=inv([A{j}*V{j} Pi]);
+    
+    tmp=norm(X*B{j} - [F1{j}; F2{j}]);
+    if tmp > diff
+        diff=tmp;
+    end
+    
+    tmp=norm(X*Psi{j} - [G1{j}; G2{j}]);
+    if (tmp > diff)
+        diff=tmp;
+    end
+end
+
+for i=1:h
+    X=zeros(s,n);
+    for j=1:h
+        X=X+P(i,j)*F2{j};
+    end
+    tmp=norm(X*V{i});
+    if tmp > diff
+        diff=tmp;
+    end;
+end