-
Notifications
You must be signed in to change notification settings - Fork 0
/
plotspglarma.m
104 lines (90 loc) · 3.36 KB
/
plotspglarma.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
function [] = plotspglarma(Y, X, link,tests,fitted,phat,sdhat,residualtype)
% four plots are currently available:
% 1) a time series plot with observed values of the dependent variable,
% fixed effects fit (SP-GLM), and GLARMA (SP-GLARMA)fit;
% 2) an ACF plot of residuals;
% 3) a plot of residuals against time;
% 4) the PIT histogram;
%
% By default, all 4 plots are provided.
%
% fitted - fitted values of the SP-GLARMA fit
% phat - phat output from SP-GLARMA
% sdhat - sdhat output form SP-GLARMA
% residualtype - type of the residual - Pearson / Score (for getting the
% right title)
switch residualtype
case 'Pearson'
[~, ~, fitted_fix, ~, ~, ~] = spglm4(Y,X, link, tests); % fitted effects fit
n = length(Y);
% ts plot
subplot(2,2,1)
h = plot(1:n,Y, 1:n,fitted_fix,1:n,fitted);
title('Observed vs Fixed vs GLARMA','FontSize',10);
set(h(1),'LineWidth',1.5,'color','k','LineStyle', '-');
set(h(2),'LineWidth',1.5,'color','red');
set(h(2),'LineWidth',1.5,'color','blue');
ymin = min([Y;fitted_fix;fitted]);
ymax = max([Y;fitted_fix;fitted]);
axis([1 n ymin ymax]);
ylabel('Y','FontSize',10)
xlabel('Time','FontSize',10)
legend('Observed','Fixed','GLARMA')
% acf plot
subplot(2,2,2)
resid = (Y-fitted)./sdhat;
autocorr(resid)
title('ACF of Pearson Residuals','FontSize',10);
ylabel('Sample ACF','FontSize',10);
xlabel('Lag','FontSize',10);
% resid plot
subplot(2,2,3)
plot(1:n,resid)
axis([1 n min(resid) max(resid)]);
title('Pearson Residuals','FontSize',10);
ylabel('Residuals','FontSize',10)
xlabel('Time','FontSize',10)
% PIT plot
subplot(2,2,4)
SPPIT_hist(Y,phat);
title('PIT for SP GLARMA','FontSize',10);
xlabel('Probability Integral Transform','FontSize',10)
ylabel('Relative Frequency','FontSize',10)
case 'Score'
[~, ~, fitted_fix, ~, ~, ~] = spglm4(Y,X,link, tests); % fitted effects fit
n = length(Y);
% ts plot
subplot(2,2,1)
h = plot(1:n,Y, 1:n,fitted_fix,1:n,fitted);
title('Observed vs Fixed vs GLARMA','FontSize',10);
set(h(1),'LineWidth',1.5,'color','k','LineStyle', '-');
set(h(2),'LineWidth',1.5,'color','red');
set(h(2),'LineWidth',1.5,'color','blue');
ymin = min([Y;fitted_fix;fitted]);
ymax = max([Y;fitted_fix;fitted]);
axis([1 n ymin ymax]);
ylabel('Y','FontSize',10)
xlabel('Time','FontSize',10)
% acf plot
subplot(2,2,2)
resid = (Y-fitted)./sdhat;
autocorr(resid)
title('ACF of Score Residuals','FontSize',10);
ylabel('Sample ACF','FontSize',10);
xlabel('Lag','FontSize',10);
% resid plot
subplot(2,2,3)
plot(1:n,resid)
axis([1 n min(resid) max(resid)]);
title('Score Residuals','FontSize',10);
ylabel('Residuals','FontSize',10)
xlabel('Time','FontSize',10)
% PIT plot
subplot(2,2,4)
SPPIT_hist(Y,phat);
title('PIT for SP GLARMA','FontSize',10);
xlabel('Probability Integral Transform','FontSize',10)
ylabel('Relative Frequency','FontSize',10)
otherwise
warning('Unexpected residual type')
end