Bug fixed

Fixed a bug in prediction. The solver does not seem to run well.

Author: junfengwen

RarmaFcns.m

@@ -42,10 +42,19 @@
         % Does not use Z^(j)_{:, T-madim+j+1:T} 
         result = zeros(xdim, numsamples);
         T = size(B,2);
+%         % JF: not working. This is a bug I never noticed, because at a
+%         time, we decided that the prediction (in
+%         convex_multi_view_model.m, line 300) should not use the "noise"
+%         for future prediction.
+%         for j = 1:madim
+%             idx = RarmaFcns.blockInd(j,xdim);
+%             Tidx = T:-1:max(T-numsamples+1,j);
+%             result(:,Tidx) = result(:,Tidx) + B(idx, Tidx-j+1);           
+%         end
         for j = 1:madim
             idx = RarmaFcns.blockInd(j,xdim);
-            Tidx = T:-1:max(T-numsamples+1,j);
-            result(:,Tidx) = result(:,Tidx) + B(idx, Tidx-j+1);           
+            Tidx = numsamples:-1:max(1,numsamples-T+j);
+            result(:,Tidx) = result(:,Tidx) + B(idx, Tidx-numsamples+T-j+1);
         end
     else
         result = zeros(xdim, numsamples);
@@ -94,7 +103,6 @@
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%% Prediction Functions %%%%%%%%%%%%%%%%%%%%%%%%%%
-
 % Model contains 
 % model.Aall, model.B, model.z and model.zparam for the exponential scalar variable
 % Phistart = [Phi_1, ..., Phi_startnum]
@@ -107,7 +115,7 @@
 % Epsilonstart can either be Zstart or the epsilon; this is
 % determined by the size of the matrix
 
-    if isempty(Phistart) 
+    if isempty(Epsilonstart) 
         Xiterated = iterate_predict_ar(Xstart, model, horizon, opts);
         return;
     end
@@ -120,12 +128,12 @@
     else  
         Xminusone = Xstart(:, (end-r+1):end);
     end  
-    Zq = Phistart(:, (end-opts.madim+1):end);
+    Zq = Epsilonstart(:, (end-opts.madim+1):end);
     if size(Zq, 1) ~= opts.madim*xdim
         if isempty(model.B)
-            Zq = zeros(opts.madim*xdim, size(Xstart,2))
+            Zq = zeros(opts.madim*xdim, size(Xstart,2));
         else    
-            Zq = model.B*Phistart(:, (end-opts.madim+1):end);
+            Zq = model.B*Epsilonstart(:, (end-opts.madim+1):end);
         end
     end
 
@@ -147,6 +155,7 @@
 % If Phistart is empty, then compute it before proceeding
 % The predicted phi will use a generative Laplace model
 
+    xdim = size(Xstart, 1);
     Xminusone = Xstart(:, (end-opts.ardim+1):end);
     xt = RarmaFcns.computeAR(model.A, Xminusone, opts.ardim, xdim, 1);
     Xiterated = xt;

RarmaUtilities.m

@@ -31,7 +31,27 @@
                     end
                 end
             end
-
         end
+        
+        function [isStable, eigs] = checkStable(A, threshold)
+        % This function check whether a model is stable based on the
+        % eigenvalues of A
+            if nargin < 2
+                threshold = 1.0;
+            end
+            [xdim, tmp] = size(A);
+            ardim = round(tmp/xdim);
+            bigA = A;
+            if ardim < 1
+                error('ardim < 1');
+            elseif ardim ~= 1
+                bigA = [bigA; eye((ardim-1)*xdim),...
+                    zeros((ardim-1)*xdim,xdim)];
+            end
+            eigs = abs(eig(bigA));
+            maxEig = max(eigs); % should be <=1
+            isStable = (maxEig<=threshold);
+        end
+        
     end
 end

demo.m

@@ -1,8 +1,11 @@
 clear
 clc
 
-num_reps = 1;
+rng(100);
+num_reps = 2;
 Models = cell(num_reps,1);
+Xpredictall = cell(num_reps,1);
+Err = zeros(num_reps,1);
 
 %% Generate data
 opts = [];
@@ -19,7 +22,12 @@
 opts.reg_wgt_ma = 1e-1;
 for ii = 1:num_reps
     Models{ii} = rarma(Xtrainall{ii},opts);
+    [isStable, eigs] = RarmaUtilities.checkStable(Models{ii}.A)
 end
 
 %% Prediction and Evaluation
-
+for ii = 1:num_reps
+    Xpredictall{ii} = Models{ii}.predict(Xtrainall{ii},...
+        size(Xtestall{ii},2), opts);
+    Err(ii) = sum(sum((Xpredictall{ii}-Xtestall{ii}).^2));
+end

rarma.m

@@ -126,7 +126,8 @@
   if opts.recover == 1
       [model.B, model.Epsilon] = recoverModels(Z);
   end
-  model.predict = @(Xstart, Phistart, horizon, opts)(RarmaFcns.iterate_predict(Xstart, model, horizon, opts));    
+%   model.predict = @(Xstart, Epsilonstart, horizon, opts)(RarmaFcns.iterate_predict(Xstart, Epsilonstart, model, horizon, opts));
+  model.predict = @(Xstart, horizon, opts)(RarmaFcns.iterate_predict_ar(Xstart, model, horizon, opts));
 end
 
 function [A, f] = iterateA(A, Z)