%% Intercept Identifiability Illustration
%% Specify the regARIMA model

Mdl0 = regARIMA('D',1,'AR',{0.8 -0.4},'MA',0.3,...
    'Intercept',0.5,'Variance', 0.2);

%
%% Simulate 100 observations

rng(1);
T = 1000;            
y = simulate(Mdl0, T);

%
%% Fit the model

Mdl = regARIMA('ARLags',1:2,'MALags',1,'D',1);...
    % "Empty" model to pass into estimate
[EstMdl,EstParamCov] = estimate(Mdl,y,'Display','params');

%
%% Plot the intercept profile likelihood

c = linspace(Mdl0.Intercept - 50,...
    Mdl0.Intercept + 50,100); % Grid of intercepts
logL = nan(numel(c),1); % For preallocation

for i = 1:numel(logL)
    EstMdl.Intercept = c(i);
    [~,~,~,logL(i)] = infer(EstMdl,y);
end

figure
plot(c,logL)
title('Profile Log-Likelihood with Respect to the Intercept')
xlabel('Intercept')
ylabel('Loglikelihood')

% Copyright 2015 The MathWorks, Inc.