Forecasting function

Give the following function a try. Let me know if it works.

function YP = forecast(model,data,K, Init)
%FORECAST Forecast a time series K steps into the future
%
% YF = FORECAST(MODEL,DATA,K)
%
% DATA: Existing data up to time N, an IDDATA object.
% MODEL: The model as any IDMODEL object, IDPOLY, IDSS, IDARX or IDGREY.
% K: The time horizon of forecasting, a positive integer with the number of 
samples
% YF: The forecasted output after time N, an IDDATA object with output
% only, covering the time span N+1:N+K.
%
%  YF = FORECAST(MODEL,DATA,K, INIT)
%  wehere INIT is 'z' or 'e' allows specification of initial conditions (at
%  time = Data.SamplingInstants(1)).
%
% See also idmodel/predict, which computes a fixed horizon prediction
% along the existing data record.

[N, ny] = size(data); % assume data is iddata
Mss = idss(model);
ord = size(pvget(Mss,'A'),1);
if ord>N
 error('Forecast:TooFewSamples','The data should contain at least %d 
samples.',ord)
end

if nargin<4, Init = 'e'; end

yp = zeros(K,ny);
mp = getPredictor(Mss);
[Ap,Bp,Cp] = ssdata(mp);
if Init=='z'
 xt = ltitr(Ap, Bp, data.y); % use zero init
 x0 = xt(end,:);
else % Init == 'e'
 [A1,B1,C1,D1,K1] = ssdata(Mss);
 x00 = 
x0est(data.y,A1,B1,C1,D1,K1,size(C1,1),size(B1,2),250e3,eye(size(C1,1)));
 x0 = ltitr(Ap,Bp,data.y,x00); x0 = x0(end,:);
end

u = [data.y(end,:); zeros(1,ny)];
for ct = 1:K
 xt = ltitr(Ap, Bp, u, x0);
 x0 = xt(end,:);
 yp(ct,:) = (Cp*x0.').';
 u = [yp(ct,:); zeros(1,ny)];
end

YP = data; YP.y = yp; YP.u = []; YP.Name = '';
YP.UserData = []; YP.Notes = '';
YP.Tstart = data.Ts*(N+1);

%------local function -----------------------------------
function mp = getPredictor(sysd)

[A,B,C,D,K] = ssdata(sysd);
Ts = sysd.Ts;
Ny = size(D,1);
mp = idss(A-K*C, [K B-K*D], C, [zeros(Ny), D], zeros(size(A,1),Ny), 'Ts', 
Ts);
mp.InputDelay = [zeros(Ny,1); sysd.InputDelay];

Dear Rajiv, thanks for sharing that function. I tried and for the sine there was no problem I used the same script to forecast the value

x = 0:0.1:30;
y = sin(x);
plot(y,'ob')
ndata = length(y);
y = y';
ts = 1;
data = iddata(y,[],ts);
na = 2;

model = ar(data,na);
K = 1;
cont = 1;
limite = 100;
while(cont<limite)
YP = forecastRaij(model,data,K,'z');
y = [data.y;YP.y];
data = iddata(y,[],ts);
cont = cont + 1;
end

YP = forecastRaij(model,data,K,'z');
y = [data.y;YP.y];
data = iddata(y,[],ts);

hold on
plot(y,'r')

%Code that verify the values forecasted
x = 0:0.1:((ndata+limite-1)/10);
y = sin(x);

y=y';

for ii = 1:((ndata+limite-1)/10)
    err(ii,1) = y(ii,1)-data.y(ii,1);
end
figure(2)
plot(err)

Let me know if there is somenthing that is not right. But When I tried with my data the value are not right, I upload the .dat file that I'm using here http://www.mediafire.com/?5yhpruipa8qp849 . The code that I use is the following:

Day = load('MyDataCorMean.dat');
% y = Day(:,1); %solar
days= 4;
y = Day((1:288*days),1); %solar
ndata = length(y);
ts = 1;
data = iddata(y,[],ts);
na = 2;
nc = 2;

model = ar(data,na);
K = 1;
cont = 1;
limit = 288*4;
YP = predict(model,data,K);
plot(YP.y,'r')
hold on
plot(y)
while(cont<limit)
YP = forecastRaij(model,data,K,'z');
y = [data.y;YP.y];
data = iddata(y,[],ts);
cont = cont + 1;
end

plot(data.y)

Mario, You are not using forecast correctly. The idea is to NOT call it in a loop (that defeats the purpose). You should simply do:

YP = forecastRaij(model,data,limit,'z');

The large prediction error (see model.es.FPE) means that the model is probably not good. One step ahead prediction is not a big task even for a poor model (just try the model y(t) = y(t-1) which wouldn't do too bad either).

Finally note that your model is a stable autoregressive process with fading memory. Any forecasting using such a model would eventually (rather quickly for your model) go to zero values. If you start with zero initial values then forecasted values would also be zero for an AR process because it has finite memory (two samples in your case). You might need an ARMA model (rather than an AR model) and perhaps a persisting process (root of "A" polynomial at unity) as in the case of seasonal ARMA models.

I would recommend investing some time into understanding what prediction and forecasting are about so that you can make educated modeling assumptions.

Rajiv thank you for your continue advices I really appreciate. I understood why I used in a not proper way your function. Now the problem is that maybe is not possible to realize a forecast for my data. I made a for loop to evaluate the value for the ARMA model but the best result is an error of 2.xxx instead the error of the model of the sine is somenthing like xxx*e-30.

na = 10;
nc = 10;
 er = zeros(na,nc);
 for ii=1:na
     for jj=1:nc
         model = armax(data,[ii jj]);
         er(ii,jj) = model.es.FPE;
     end
 end

My only problem is to find a good model and after that I can use your function and forecast the data that I want. But how can I found such a model? Thanks

In R2012a, you can also use a new "FORECAST" function. See help on "idParametric/forecast".