001 function [tpe,Pout,I,tpe0] = tpextrapolate(tp,N,Pin,plotflag)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
052
053
054
055
056
057
058
059
060
061
062
063
064
065 ni = nargin;
066 no = nargout;
067 error(nargchk(2,4,ni));
068
069 if ni<3, Pin=[]; end
070 if ni<4, plotflag=[]; end
071
072 if isempty(plotflag)
073 plotflag=0;
074 end
075
076
077
078
079
080 Pout.method = 'exp';
081 Pout.LCfrac = [];
082 Pout.u_lev = [];
083 Pout.lim = [];
084
085
086 if ~isempty(Pin)
087 Fname = fieldnames(Pin);
088 for i = 1:length(Fname)
089 Pout = setfield(Pout,Fname{i},getfield(Pin,Fname{i}));
090 end
091 end
092
093 method = lower(Pout.method);
094 if ~(strcmp(method,'exp') | strcmp(method,'gpd'))
095 error(['Undefined method: ' Pout.method]),
096 end
097
098 if isempty(Pout.u_lev)
099 lc = tp2lc(tp);
100
101 LCmax=max(lc(:,2));
102 if isempty(Pout.LCfrac)
103 Pout.LCfrac = 1/sqrt(LCmax);
104 end
105
106 nLCfrac = LCmax*Pout.LCfrac;
107 imax = max(find(lc(:,2)>nLCfrac));
108 imin = min(find(lc(:,2)>nLCfrac));
109 umax = lc(imax,1);
110 umin = lc(imin,1);
111 Pout.u_lev = [umin umax];
112 end
113
114 u_min = Pout.u_lev(1);
115 u_max = Pout.u_lev(2);
116
117 if tp(1,end)>tp(2,end)
118 StartMax=1; StartMin=2;
119 else
120 StartMax=2; StartMin=1;
121 end
122
123
124
125 Imax = 2*(find(tp(StartMax:2:end,end)>u_max)-1)+StartMax;
126 Imin = 2*(find(tp(StartMin:2:end,end)<u_min)-1)+StartMin;
127
128 y_max = tp(Imax,end)-u_max;
129 y_min = -(tp(Imin,end)-u_min);
130
131 a_max = mean(y_max);
132 a_min = mean(y_min);
133 if strcmp(method,'gpd')
134 [gpd_max,cov_max] = wgpdfit(y_max,'ml',0);
135 [gpd_min,cov_min] = wgpdfit(y_min,'ml',0);
136 ci_k_max = [gpd_max(1)-2*sqrt(cov_max(1)) gpd_max(1)+2*sqrt(cov_max(1))]
137 ci_k_min = [gpd_min(1)-2*sqrt(cov_min(1)) gpd_min(1)+2*sqrt(cov_min(1))]
138 end
139
140 if plotflag>0
141 subplot(2,2,1),plot(1:length(y_max),y_max,'.'),
142 title(['Exceedances of maxima above u_{max}=' num2str(u_max)]), ylabel('Exceedances of maxima')
143 subplot(2,2,2), if length(y_max)>0, wexpplot(y_max), end
144 subplot(2,2,3),plot(1:length(y_min),y_min,'.'),
145 title(['Exceedances of minima below u_{min}=' num2str(u_min)]), ylabel('Exceedances of minima')
146 subplot(2,2,4), if length(y_min)>0, wexpplot(y_min), end
147 drawnow
148
149 Pout.Zmin = y_min;
150 Pout.Zmax = y_max;
151 end
152
153 if N>0
154 meth=2;
155 conservarive = 1;
156
157 if meth==2
158 [yy_max,IImax] = sort(y_max);
159 [yy_min,IImin] = sort(y_min);
160 end
161 tpe = [];
162 for k = 1:N
163
164
165
166 if length(y_max)>0
167 if strcmp(method,'exp')
168 yr_max = wexprnd(a_max,length(y_max),1);
169 else
170 yr_max = wgpdrnd(gpd_max(1),gpd_max(2),0,length(y_max),1);
171 end
172 else
173 yr_max=[];
174 end
175
176 if length(y_min)>0
177 if strcmp(method,'exp')
178 yr_min = wexprnd(a_min,length(y_min),1);
179 else
180 yr_min = wgpdrnd(gpd_min(1),gpd_min(2),0,length(y_min),1);
181 end
182 else
183 yr_min=[];
184 end
185
186 if meth ==1
187
188
189 tpr = tp;
190 tpr(Imax,end) = u_max + yr_max;
191 tpr(Imin,end) = u_min - yr_min;
192 else
193
194
195
196
197
198 yr_max = sort(yr_max);
199 yr_min = sort(yr_min);
200
201 tpr = tp;
202 tpr(Imax(IImax),end) = u_max+yr_max;
203 tpr(Imin(IImin),end) = u_min-yr_min;
204 end
205
206 if conservarive
207 I=find(tpr(Imin,end)>tp(Imin,end));
208 tpr(Imin(I),end) = tp(Imin(I),end);
209 I=find(tpr(Imax,end)<tp(Imax,end));
210 tpr(Imax(I),end) = tp(Imax(I),end);
211 end
212
213 tpe = [tpe; tpr];
214 end
215
216 if no>3
217 tpe0=tpe;
218 end
219
220
221 if ~isempty(Pout.lim)
222 minlim=Pout.lim(1); maxlim=Pout.lim(2);
223 if isnumeric(maxlim)
224 I = find(tpe(:,end)>maxlim);
225 if ~isempty(I)
226 tpe(I,end)=maxlim;
227 end
228 end
229 if isnumeric(minlim)
230 I = find(tpe(:,end)<minlim);
231 if ~isempty(I)
232 tpe(I,end)=minlim;
233 end
234 end
235 end
236
237
238 tpe = rfcfilter(tpe,0,1);
239
240
241 else
242 tpe = [];
243 end
244
245
246 Pout.a_min = a_min;
247 Pout.a_max = a_max;
248 Pout.n_min = length(y_min);
249 Pout.n_max = length(y_max);
250
251 if strcmp(method,'gpd')
252 Pout.gpd_min = gpd_min;
253 Pout.gpd_max = gpd_max;
254 Pout.ci_k_min = ci_k_min;
255 Pout.ci_k_max = ci_k_max;
256 end
257
258 if N>0
259
260 if no>2
261 if tpe(1,end)>tpe(2,end)
262 StartMax=1; StartMin=2;
263 else
264 StartMax=2; StartMin=1;
265 end
266
267
268 I=[];
269 I.max = 2*(find(tpe(StartMax:2:end,end)>u_max)-1)+StartMax;
270 I.min = 2*(find(tpe(StartMin:2:end,end)<u_min)-1)+StartMin;
271 end
272 end
273
274
275 if plotflag>1
276
277 m = mean(tp(:,end));
278
279 n=500;
280 Umin = linspace(min(tp(:,end)),m,n);
281 Umax = linspace(m,max(tp(:,end)),n);
282 Amin=zeros(1,n); Amax=zeros(1,n);
283 dAmin=zeros(1,n); dAmax=zeros(1,n);
284 for k=1:n
285 Imax = 2*(find(tp(StartMax:2:end,end)>Umax(k))-1)+StartMax;
286 Imin = 2*(find(tp(StartMin:2:end,end)<Umin(k))-1)+StartMin;
287
288 y_max = tp(Imax,end)-Umax(k);
289 if ~isempty(y_max), Amax(k) = mean(y_max); dAmax(k)=Amax(k)/sqrt(length(y_max)); end
290 y_min = -(tp(Imin,end)-Umin(k));
291 if ~isempty(y_min), Amin(k) = mean(y_min); dAmin(k)=Amin(k)/sqrt(length(y_min)); end
292 end
293
294 figure
295 subplot(2,1,1)
296 plot(Umax,Amax,'r',Umax,Amax-2*dAmax,'b:',Umax,Amax+2*dAmax,'b:'), grid
297 title('Maxima - Choose a level where the estimate is stable'),
298 ylabel('Estimated mean exceedance, m'), xlabel('Upper threshold level, u_{max}')
299 subplot(2,1,2)
300 plot(-Umin,Amin,'r',-Umin,Amin-2*dAmin,'b:',-Umin,Amin+2*dAmin,'b:'), grid
301 title('Minima - Choose a level where the estimate is stable'),
302 ylabel('Estimated mean exceedance, m'), xlabel('Lower threshold level, -u_{min}')
303
304 Pout.Umin = Umin;
305 Pout.Amin = Amin;
306 Pout.Umax = Umax;
307 Pout.Amax = Amax;
308 end
309
310