%% Estimate Linear Model Coefficients Using Iterative FGLS
% This example expands on the analysis in Estimate FGLS Coefficients of Models Containing ARMA Errors. Create this regression model with ARMA(1,2) errors, where \varepsilon_t is Gaussian with mean 0 and variance 1.
% \begin{array}{*{20}{c}}
% {{y_t} = 1 + {x_t}\left[ {\begin{array}{*{20}{c}}
% 2\\
% 3
% \end{array}} \right] + {u_t}}\\
% {{u_t} = 0.6{u_{t - 1}} + {\varepsilon _t}} - 0.3\varepsilon_{t-1} +
% 0.1\varepsilon_{t-1}.
% \end{array}

beta = [2 3];
phi = 0.2;
theta = [-0.3 0.1];
Mdl = regARIMA('AR',phi,'MA',theta,'Intercept',1,'Beta',beta,'Variance',1);

%
% Simulate 500 periods of 2-D standard Gaussian values for x_t, and then simulate responses using Mdl.

numObs = 500;
rng(1); % For reproducibility
X = randn(numObs,2);
y = simulate(Mdl,numObs,'X',X);

%
% Convert the ARMA model polynomial to an infinite-lag AR model polynomial using arma2ar. By default, arma2ar returns the coefficients for the first 10 terms. Find the number of terms that exceed 0.00001.

arParams = arma2ar(phi,theta);
arLags = sum(abs(arParams) > 0.00001);

%
% Estimate the regression coefficients using three iterations of FGLS, and specify the number of lags in the AR innovation model (arLags). Also, specify to plot the coefficient estimates and their standard errors for each iteration, and to display the final estimates and the OLS estimates in tabular form.

[coeff,~,EstCoeffCov] = fgls(X,y,'innovMdl','AR','arLags',arLags,...
    'numIter',3,'plot',{'coeff','se'},'display','final');

% The algorithm seems to converge after the first iteration, and the estimates are close to the OLS estimates, with the standard errors being slightly smaller.
% Properties of iterative FGLS estimates in finite samples are difficult to establish. For asymptotic properties, one iteration of FGLS is sufficient. fgls supports iterative FGLS for experimentation.
% If the estimates or standard errors show instability after successive iterations, then the estimated innovations covariance might be ill conditioned. Consider scaling the residuals using the 'resCond' name-value pair argument to improve the conditioning of the estimated innovations covariance.
% Copyright 2015 The MathWorks, Inc.