Home > wafo > wdemos > rfcdemo1.m

# rfcdemo1

## PURPOSE

Demo for switching AR(1)-processes.

## SYNOPSIS

[F_RFC] = refdemo1(demoNr,P,A,m,s2,param)

## DESCRIPTION

```  RFCDEMO1  Demo for switching AR(1)-processes.

[F_RFC] = rfcdemo1;
[F_RFC] = rfcdemo1(demoNr);
[F_RFC] = rfcdemo1(demoNr,P);
[F_RFC] = rfcdemo1(demoNr,P,A,m,s2,param);

Input:
demoNr = 1 / 2 / 3 / 0
Default: 1. 0 = Define your own process.
P      = Transition matrix for regime process.  [rxr]
A      = AR coefficients.                       [rx2]
m      = Mean koefficients.                     [rx1]
s2     = Innovation variances.                  [rx1]
param  = Parameter vector, [a b n], defines discretization.

Output:
F_RFC  = Calculated rainflow intensity          [nxn]

Demo for switching AR(1)-processes, by using
switching Markov chains.
See Examples 3.1-3.3 in PhD-thesis:
P. Johannesson (1999): Rainflow Analysis of
[http://www.maths.lth.se/matstat/staff/pj/PhD/]

Example:
rfcdemo1;
rfcdemo1(2);
P=[0.98 0.02; 0.05 0.95]; rfcdemo1(1,P);

## CROSS-REFERENCE INFORMATION

This function calls:
 armaspec Computes the spectral density for an AR- or ARMA-model. cc2lc Calculates the number of upcrossings from a cycle count cmatplot Plots a cycle matrix, e.g. a rainflow matrix. cocc Plots cycles as points together with isolines of a cycle matrix. dat2tp Extracts turning points from data, discar Discretizes an AR(1) process. hmmplot plots a Hidden Markov Model. levels Calculates discrete levels given the parameter matrix. mc2rfm Calculates the rainflow matrix/intensity for a Markov chain. mc2stat Calculates the stationary distribution for a Markov chain. sarmasim Simulates a switching ARMA-process. smc2rfm Calculates the rainflow matrix/intensity for a switching Markov chain. tp2rfc Finds the rainflow cycles from the sequence of turning points. waforoot Root directory of WAFO installation. addpath Add directory to search path. axis Control axis scaling and appearance. delete Delete file or graphics object. drawnow Flush pending graphics events. error Display message and abort function. exist Check if variables or functions are defined. figure Create figure window. filesep Directory separator for this platform. gcf Get handle to current figure. grid Grid lines. hold Hold current graph. num2str Convert number to string. (Fast version) pause Wait for user response. plot Linear plot. rmpath Remove directory from search path. set Set object properties. square Square wave generation. stairs Stairstep plot. subplot Create axes in tiled positions. title Graph title. warning Display warning message; disable or enable warning messages. xlabel X-axis label. ylabel Y-axis label.
This function is called by:

## SOURCE CODE

```001 function [F_RFC] = refdemo1(demoNr,P,A,m,s2,param)
002
003 % RFCDEMO1  Demo for switching AR(1)-processes.
004 %
005 % [F_RFC] = rfcdemo1;
006 % [F_RFC] = rfcdemo1(demoNr);
007 % [F_RFC] = rfcdemo1(demoNr,P);
008 % [F_RFC] = rfcdemo1(demoNr,P,A,m,s2,param);
009 %
010 % Input:
011 % demoNr = 1 / 2 / 3 / 0
012 %          Default: 1. 0 = Define your own process.
013 % P      = Transition matrix for regime process.  [rxr]
014 % A      = AR coefficients.                       [rx2]
015 % m      = Mean koefficients.                     [rx1]
016 % s2     = Innovation variances.                  [rx1]
017 % param  = Parameter vector, [a b n], defines discretization.
018 %
019 % Output:
020 % F_RFC  = Calculated rainflow intensity          [nxn]
021 %
022 %   Demo for switching AR(1)-processes, by using
023 %   switching Markov chains.
024 %   See Examples 3.1-3.3 in PhD-thesis:
025 %     P. Johannesson (1999): Rainflow Analysis of
027 %     [http://www.maths.lth.se/matstat/staff/pj/PhD/]
028 %
029 % Example:
030 %   rfcdemo1;
031 %   rfcdemo1(2);
032 %   P=[0.98 0.02; 0.05 0.95]; rfcdemo1(1,P);
033 %
035
036 % Tested  on Matlab  5.3
037 %
038 % History:
039 % Created by PJ (Pär Johannesson) 1997
040 % Updated by PJ 19-May-2000
041 %   updated for WAFO
042 % Updated by PJ 07-Jul-2005
043 %   warning off
044
045 warning_state=warning;  % Get current warning-settings
046 warning off
047
048 fprintf(1,'Welcome to demo 1 for "Rainflow Cycles for Switching Processes"!\n');
049 fprintf(1,'It demonstrates calculation of the rainflow intensity\n');
050 fprintf(1,'for a switching AR(1)-process.\n');
051 fprintf(1,'Type "help rfcdemo1" for further information\n\n');
052
053 % Check input arguments
054
055 ni = nargin;
056 no = nargout;
057 error(nargchk(0,6,ni));
058
059 Zstr='123456789';
060
061 % Add path to rfcdemo1
062 demoDir = [waforoot filesep 'wdemos' filesep 'rfcdemo1'];
064
065 % Delete all figure windows
066
067 slut=0;
068 h=gcf;
069 while ~slut
070   h_prev=h;
071   delete(h);
072   h=gcf;
073   if h==h_prev, slut=1; end
074 end
075
076
077 if ni == 0, demoNr=1; end
078
079 % Define switching AR(1)-process
080
081 if demoNr == 1
082   p1=0.02; p2=0.01;
083   P0 = [1-p1 p1; p2 1-p2];
084   A = [1 -0.5; 1 0.5];
085   m = 2*[-0.5 1.5]';
086   s2 = [1 1]';
087   param = [-8 8 64];
088 elseif demoNr == 2
089   p1=0.01; p2=0.02;
090   P0 = [1-p1 p1; p2 1-p2];
091   A = [1 -0.5; 1 0.5];
092   m = [0 0]';
093   s2 = [4 1]';
094   param = [-10 10 64];
095 elseif demoNr == 3
096   p12=0.02; p13=0.005;
097   p21=0.01; p23=0.01;
098   p31=0.005; p32=0.02;
099   P0 = [1-p12-p13 p12 p13; p21 1-p21-p23 p23; p31 p32 1-p31-p32];
100   A = [1 -0.5; 1 -0.3; 1 0.5];
101   m = 2*[-0.5 0 1.5]';
102   s2 = [1 1 1.44]';
103   param = [-8 8 64];
104 elseif demoNr == 0
105   if ni < 6
106     error('demoNr=0: Too few input arguments.');
107   end
108 end
109
110 if ~exist('P')
111   P=P0;
112 end
113
114 % Initiation
115
116 u=levels(param);    % Discretization levels
117 delta = u(2)-u(1);  % Discretization step
118 r=length(P);        % Number of regime states
119 n=param(3);         % Number of discretization levels
120 C=ones(r,1);        %
121
122 % Calculate statistics: mean , std, proportion, mean time
123
124 mX = m./sum(A')';
125 sX = sqrt(s2./(1-A(:,2).^2));
126 statP = mc2stat(P);
127 tau = 1./(1-diag(P))';
128
129 % Parameters
130
131 fprintf('Parameters\n');
132 fprintf(' z      a(z)      m(z)      s(z)\n');
133 fprintf('%2d  %8g  %8g  %8g\n',[(1:r)' A(:,2) m sqrt(s2)]')
134 P
135 fprintf('\n')
136
137 % Statistics
138 fprintf('Statistics\n');
139 fprintf('\n z    m_X(z)    s_X(z)     ro(z)    tau(z)\n');
140 fprintf('%2d  %8g  %8g  %8g  %8g\n',[(1:r)' mX sX statP' tau']')
141 fprintf('\n')
142
143 % Calculate spectral densities
144
145 fprintf(1,' -- Calculate spectral densities.\n');
146
147 set(gcf,'Name','Power spectra')
148
149 for i=1:r
150   Si = armaspec(1,A(i,:),s2(i,:),100);
151   subplot(r,1,i)
152   plot(Si(:,1),Si(:,2)), grid
153   ylabel(['S(f), regime ' Zstr(i)])
154 end
155 xlabel('frequency, f')
156
157 drawnow, disp('Hit any key to continue.'); pause;
158
159 % Simulate MARX(2,1)-process
160
161 fprintf(1,' -- Simulating.\n');
162
163 T=10000;   % Length of simulation
164 Tinit=500; % Length of initiation (to remove transients)
165
166 [x,z,e] = sarmasim(C,A,m,s2,P,T,Tinit);
167
168 figure
169 set(gcf,'Name','Switching AR(1)-process')
170
171 I = 1:501;
172 hmmplot(x(I),z(I),I-1,[1 r],'','Switching Process, X(k)',1);
173 ylabel('regime, Z(k)'),xlabel('time, k')
174
175 drawnow, disp('Hit any key to continue.'); pause;
176
177 %
178 % Discretize AR(1)-processes
179 % Approximate with Markov chain
180 %
181
182 ud = u; ud(1)=u(1)-0.25*(u(n)-u(1)); ud(n)=u(n)+0.25*(u(n)-u(1));
183
184 figure
185 set(gcf,'Name','Markov chain for Discretized AR(1)-processes')
186
187 if demoNr == 1 | demoNr ==2 | demoNr ==3
189
191   for i=1:r
192     subplot(1,r,i)
193     eval(['Q{i}=Q' num2str(i) ';']);
194     cmatplot(uu,uu,Q{i},3), axis square
195     title(['Regime ' Zstr(i)])
196   end;
197   drawnow
198 else
199
200   fprintf(1,' -- Discretizing AR-processes. (This can take some time!)\n');
201   for i=1:r
202     fprintf(1,' -- Discretizing process %d.\n',i);
203     a1i = A(i,2);
204     si = sqrt(s2(i));
205     mi = mX(i);
206     [Qi,uu] = discar(ud,a1i,mi,si);
207     Q{i}=Qi;
208     %eval(['Q{i}' num2str(i) '=Qi;']);
209     subplot(r,1,i)
210     cmatplot(uu,uu,Qi,3), axis square
211     title(['Regime ' Zstr(i)])
212     drawnow
213   end
214 end
215
216 drawnow, disp('Hit any key to continue.'); pause;
217
218 %
219 % Compute RFC counting distribution for SMC
220 %
221
222 fprintf(1,' -- Calculating rainflow intensity.\n');
223
224 Qstr='';
225 for i=1:r
226   Qstr = [Qstr ',Q' Zstr(i)];
227 end
228
229 [F_RFC,mu_RFC] = smc2rfm(P,Q,[2 delta]);    % definition 2
230 %eval( ['[F_RFC,mu_RFC] = smc2rfc(P,[2 delta]' Qstr ');'] );    % definition 2
231
232 for i=1:r
233   %eval( ['[F_RFC' Zstr(i) ',mu_RFC' Zstr(i) '] = mc2rfc(Q' Zstr(i) ',[2 delta]);'] );     % Regime i
234   [F_RFC0{i},mu_RFC0{i}] = mc2rfm(Q{i},[2 delta]);     % Regime i
235 end
236
237 figure
238
239 cmatplot(u,u,F_RFC,2);
240 set(gcf,'Name','Rainflow intensity for switching AR(1)-process')
241 title('3D-plot')
242
243 drawnow, disp('Hit any key to continue.'); pause;
244
245 figure
246
247 % Calculate Cycles and crossings
248
249 fprintf(1,' -- Cycles and crossings for simulated load.\n');
250
251 TP = dat2tp([(1:T)' x]);
252 RFC = tp2rfc(TP);
253 cross = cc2lc(RFC);
254
255 cocc(param,RFC,F_RFC)
256 set(gcf,'Name','Rainflow intensity for switching AR(1)-process')
257 title('Iso-lines compared with observed cycles')
258 xlabel('min'), ylabel('Max')
259
260 drawnow, disp('Hit any key to continue.'); pause;
261
262 % Plot Crossing intensity
263
264 figure
265 set(gcf,'Name','Crossing intensity')
266
267 subplot(1,2,1)
268 plot(u,diag(mu_RFC),'g--'), hold on
269 stairs(cross(:,1),cross(:,2)/T), hold off
270 title('Calculated vs. observed')
271 xlabel('level, u'); ylabel('Crossing intensity, mu(u)');
272
273 subplot(1,2,2)
274 plot(u,diag(mu_RFC),'g'), hold on
275 cross_sum=zeros(size(diag(mu_RFC)));
276 for i=1:r
277   crossi=diag(mu_RFC0{i});
278   plot(u,statP(i)*crossi,'r-.')
279   cross_sum=cross_sum+statP(i)*crossi;
280 end
281 plot(u,cross_sum,'y--'), hold off
282 title('Calculated vs. individual')
283 xlabel('level, u'); ylabel('Crossing intensity, mu(u)');
284
285 fprintf(1,'\nThank you for running this demo!\n');
286
287 % Remove path to rfcdemo1
288 demoDir = [waforoot filesep 'wdemos' filesep 'rfcdemo1'];
289 rmpath(demoDir);
290
291 warning(warning_state);
292
293```

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group

Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003