Fix cases when p=0 but q>0

Note that prediction is not meaningful when p=0

Author: junfengwen

RarmaFcns.m

@@ -30,10 +30,14 @@
 %         return;
 %     end
 
-    Xhist = RarmaFcns.generate_history(Xminusone, ardim);
-    result = A*Xhist;
-    result = [zeros(xdim, ardim) result];
-    if numsamples < size(result,2), result = result(:, end-numsamples+1:end); end
+    if ardim > 0
+        Xhist = RarmaFcns.generate_history(Xminusone, ardim);
+        result = A*Xhist;
+        result = [zeros(xdim, ardim) result];
+        if numsamples < size(result,2), result = result(:, end-numsamples+1:end); end
+    else
+        result = zeros(xdim, numsamples);
+    end
 end
 
 function result = computeMA(B, Epsilon, madim, xdim, numsamples)
@@ -42,15 +46,7 @@
         % 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 = numsamples:-1:max(1,numsamples-T+j);
@@ -115,8 +111,8 @@
 % Epsilonstart can either be Zstart or the epsilon; this is
 % determined by the size of the matrix
 
-    if isempty(Epsilonstart) 
-        Xiterated = iterate_predict_ar(Xstart, model, horizon, opts);
+    if isempty(Epsilonstart)
+        Xiterated = RarmaFcns.iterate_predict_ar(Xstart, model, horizon, opts);
         return;
     end
 

demo.m

@@ -3,17 +3,15 @@
 
 rng(10);
 num_reps = 10;
-Models = cell(num_reps,1);
-Xpredictall = cell(num_reps,1);
-isStable = zeros(num_reps,1);
-Err = zeros(num_reps,1);
 
 %% Generate data
 opts = [];
 opts.num_reps = num_reps;
 [Xstartall,Xtrainall,Xtestall] = genARMA(opts);
 
 %% Learn RARMA
+Models = cell(num_reps,1);
+isStable = zeros(num_reps,1);
 % In practice, the following paramters should be 
 % cross-validated for EACH dataset before applied
 opts = [];
@@ -23,10 +21,14 @@
 opts.reg_wgt_ma = 0.01;
 for ii = 1:num_reps
     Models{ii} = rarma(Xtrainall{ii},opts);
-    isStable(ii) = RarmaUtilities.checkStable(Models{ii}.A);
+    if opts.ardim > 0
+        isStable(ii) = RarmaUtilities.checkStable(Models{ii}.A);
+    end
 end
 
-%% Prediction and Evaluation
+%% Prediction and Evaluation (only when ardim > 0)
+Xpredictall = cell(num_reps,1);
+Err = zeros(num_reps,1);
 for ii = 1:num_reps
     Xpredictall{ii} = Models{ii}.predict(Xtrainall{ii},...
         size(Xtestall{ii},2), opts);

rarma.m

@@ -14,7 +14,7 @@
 %
 % Authors: Junfeng Wen (University of Alberta)
 %          Martha White (Indiana University) 
-%          Released: March 2015
+%          Last Update: Nov 2015
 
 if nargin < 1
   error('rarma requires at least data matrix X = [X1, ..., XT]');
@@ -46,6 +46,11 @@
   opts = RarmaUtilities.getOptions(opts, DEFAULTS);
 end
 
+% Validity check
+if opts.madim < 0 || opts.ardim < 0
+  error('Incorrect parameters');
+end
+
 % Solve for RARMA, AR or MA dependning on madim and ardim choices
 ar_solver = @solve_rarma;
 if opts.madim < 1
@@ -64,8 +69,13 @@
 % Initialize variables for learning
 [xdim, numsamples] = size(X);
 sizeA = [xdim, xdim*opts.ardim];
-Ainit = initAParams(sizeA);
 sizeZ = [xdim*opts.madim, numsamples];
+if opts.ardim > 0
+  Ainit = initAParams();
+end
+if opts.madim > 0
+  Zinit = zeros(sizeZ);
+end
 
 % START the outer optimization; save any learned variables in model along the way
 model = [];
@@ -90,32 +100,32 @@
 % Note: the code is currently not well-designed to do long-horizon prediction
 % with only a moving average models, since future moving average
 % components are not imputed; should work, however, for 1-step prediction
-  [Z, obj, iter, msg] = opts.optimizer(@(Zvec)(objZ(Zvec, zeros(sizeA))), zeros(sizeZ));
+  [Z, obj, iter, msg] = opts.optimizer(@(Zvec)(objZ(Zvec, zeros(sizeA))), Zinit(:));
   model.Z = reshape(Z, sizeZ);
   model.A = zeros(sizeA);
   if opts.recover == 1
       [model.B, model.Epsilon] = recoverModels(Z);
-  end  
-  model.predict = @(Xstart, horizon, ...
-                    opts)(RarmaFcns.iterate_predict(Xstart, [], model, horizon, opts));  
+  end
+  % Cannot really 'predict' without autoregressive component
+  model.predict = @(Xstart, horizon, opts)...
+      (RarmaFcns.iterate_predict(Xstart, [], model, horizon, opts));
 end
 
 
 function [model,obj] = solve_rarma()
 %% SOLVE_RARMA
 % Solve for A first (with Z = 0), then iterate between Z and A
-    Z = zeros(sizeZ);
-    [A, obj, iter, msg] =opts.optimizer(@(Avec)(objA(Avec, X, Z)), Ainit(:));
+    [A, obj, iter, msg] = opts.optimizer(@(Avec)(objA(Avec, X, Zinit)), Ainit(:));
     A = reshape(A, sizeA);
-    [Z, prev_obj] = iterateZ(Z,A,opts.init_stepsize);
+    [Z, prev_obj] = iterateZ(Zinit, A, opts.init_stepsize);
 
   for i = 1:opts.maxiter
     % Do A first since it returns the incorrect obj
-    A = iterateA(A,Z,opts.init_stepsize/i); % adaptive stepsize
-    [Z, obj] = iterateZ(Z,A,opts.init_stepsize/i);
+    A = iterateA(A, Z, opts.init_stepsize/i); % adaptive stepsize
+    [Z, obj] = iterateZ(Z, A, opts.init_stepsize/i);
 
     if abs(prev_obj-obj) < opts.TOL % doing minimization
-      break; 
+      break;
     end
     prev_obj = obj;
   end
@@ -189,7 +199,7 @@
   f = f + opts.reg_wgt_ma * f2;
 end
 
-function Aparams = initAParams(sizeA)
+function Aparams = initAParams()
 % INITAPARAMS
 % Could initialize in many ways; for speed, we choose
 % a regression between X = AXhist