function surx = SurrogateARMAEcon15(x,p,q,nSur)
%SURROGATEARMAECON Surrogate series using ECONOMETRICS toolbox.
%              
%	SURX=SURROGATEARMAECON(X,P,Q,NSUR) 
%   Computes NSUR surrogate series of a given time series X. It constructs 
%   the surrogates by making use of the matlab functions: GARCHFIT, 
%   GARCHSET and GARCHSIM from the ECONOMETRICS Toolbox. It first fits an 
%   ARMA(P,Q) model to series X and and then construct new samples by 
%   drawing errors from a Gaussian distribution.
%
%   INPUTS: 
%       X - a vector (the raw series).
%   Optional inputs: 
%       P,Q - integers,  orders of the ARMA process (Default: P=1,Q=0).
%       NSUR - positive integer, the number of surrogates (Default: 1).
%   OUTPUT:
%       SURX - the matrix of size (length(X) x NSUR) whose columns contain 
%              the surrogates of X.
%
%   See also: SurrogateARMABoot.

 
%   Written by L. Aguiar-Conraria and M.J. Soares
%
%   Luís AGUIAR-CONRARIA              Maria Joana SOARES                      
%   Dep. Economics                    Dep. Mathematics and Applications   
%   University of Minho               University of Minho
%   4710-057 Braga                    4710-057 Braga
%   PORTUGAL                          PORTUGAL
%                           
%   lfaguiar@eeg.uminho.pt            jsoares@math.uminho.pt 

%%%%%%%%%%%%%%%%%%%% Test for inputs and default values %%%%%%%%%%%%%%%%%%%
if ( nargin<1 )
    error('Must input series');
end

if ( nargin<2 || isempty(p) ), p=1; end
if ( nargin<3 || isempty(q) ), q=0; end
if ( nargin<4 ), nSur=1; end
    
if ( ~(isvector(x)) )
    error('x must be a vector ')
end
lx=length(x);
if (lx < p || lx< q )
	error('Length of series must be greater or equal to P and Q')
end
if ( ~(p == round(p)) || p<0 )
	 error('Model order P must be a non-negative integer')
end
if ( ~(q == round(q)) || q<0 )
	 error('Model order Q must be a non-negative integer')
end
if ( ~(nSur==round(nSur)) || nSur<=0 )
    error('nSur must be a positive integer')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mx = mean(x);       % Mean of x
sx = std(x);        % Standard deviation of x
Mx = repmat(mx,lx,nSur); % Matrix of size of surx, filled with the value mx
Sx = repmat(sx,lx,nSur); % Matrix of size of surx, filled with the value sx
x = (x - mx) ./ sx;       % Normalize x
 
% Use of functions GARCHSIM, GARCHFIT and GARCHSET to compute the
% surrogates

% [d1,d2,surx] = garchsim(garchfit(garchset('R',p,'M',q,...
%     'Display','off'),x),lx,nSur); %#ok<*ASGLU>
 
spec = arima(p,0,q); 
[ EstMdl ] = estimate(spec,x); 
surx1 = NaN(lx,nSur); 
for ii = 1 : nSur 
    surx1(:,ii) = simulate(EstMdl,lx); 
end 
 
surx = surx1;  
surx = ProcessMatrix(surx,0,1);     % Normalize marix of surrogates
surx = surx .* Sx + Mx; % Make surrogates have the same mean and standard  
                               % deviation as x
                      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end 
 
 
 
 
