Home > wafo > wsim > smcsim.m

# smcsim

## PURPOSE

Simulates a Switching Markov chain with state space.

## SYNOPSIS

[x,z] = simsmc(P,Qc,T,init);

## DESCRIPTION

` SMCSIM  Simulates a Switching Markov chain with state space.`

## CROSS-REFERENCE INFORMATION

This function calls:
 mc2stat Calculates the stationary distribution for a Markov chain. mcsim Simulates a Markov chain. error Display message and abort function. kron Kronecker tensor product. warning Display warning message; disable or enable warning messages.
This function is called by:
 test_markov Quick test of the routines in module 'markov'

## SOURCE CODE

```001 function [x,z] = simsmc(P,Qc,T,init);
002
003
004
005 %SMCSIM  Simulates a Switching Markov chain with state space.
006
007 %
008
009 % CALL: [x,z] = smcsim(P,Q,T);
010
011 %       [x,z] = smcsim(P,Q,T,init);
012
013 %
014
015 % x    = Simulated switching Markov chain
016
017 % z    = Simulated Regime process
018
019 %
020
021 % P    = Transition matrix for regime process     [rxr]
022
023 % Q      = Cell array of transition matrices        {r,1}
024
025 % Q{i}   = Transition matrix for Markov chain i     [nxn]
026
027 % T    = Length of simulation.
028
029 % init.x0 = Initial state of process x. If not given, it will start from
030
031 %          the stationary distribution of minima given z(1).
032
033 % init.z0 = Initial state of regime process. If not given, it will start
034
035 %          from the stationary distribution of the Markov chain.
036
037 %
038
039 % Simulates a Switching Markov chain with state space {1,2,...,n}.
040
041 % The regime process has state space {1,2,...,r}.
042
043 %
044
045 % Example: Simulation of a switching Markov chain with two regime states.
046
047 %   P=[0.98 0.02; 0.05 0.95]; n=16;
048
049 %   Q{1} = rand(n,n)*diag(exp(5*((n:-1:1)-1)/n));
050
051 %   Q{2} = rand(n,n)*diag(exp(5*((1:n)-1)/n)); % They will be normalized to row sum 1.
052
053 %   [x,z] = smcsim(P,Q,400); hmmplot(x,z)
054
055 %   init.z0 = 2; init.x0 = [];
056
057 %   [x,z] = smcsim(P,Q,400,init); hmmplot(x,z,[],[1 2],'','',1)
058
059 %   init.z0 = []; init.x0 = 6;
060
061 %   [x,z] = smcsim(P,Q,400,init); hmmplot(x,z,[],[1 2],'','',1)
062
063 % Example: Simulation of a Markov chain
064
065 %   P=[0.9 0.1; 0.05 0.95];
066
067 %   x = smcsim(1,P,1000);
068
069 %   plot(x)
070
071
072
073 % Tested  on Matlab  5.3
074
075 %
076
077 % History:
078
079 % Revised by PJ 19-May-2000
080
081 %   updated for WAFO
082
083 %   Corrected method for simulating starting conditions.
084
085 % Created by PJ (Pär Johannesson) 1997
086
087 %   Copyright (c) 1997 by Pär Johannesson
088
089 %   Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997
090
091
092
093 % Check input arguments
094
095
096
097 ni = nargin;
098
099 no = nargout;
100
101 error(nargchk(3,4,ni));
102
103
104
105 if ni < 4,  init = []; end
106
107
108
109 if isempty(init)
110
111   init.x0 = [];
112
113   init.z0 = [];
114
115 end
116
117
118
119 % Set constants
120
121 Zstr = '123456789';
122
123
124
125
126
127 r = length(P);     % Number of regime states
128
129 n = length(Qc{1}); % Number of levels
130
131
132
133 % Check that the rowsums of P are equal to 1
134
135
136
137 sumP = sum(P');
138
139 if sum(sumP == 1) ~= length(P)
140
141   warning(['Rowsums of P not equal to 1. Renormalizing.']);
142
143   for i = 1:length(P)
144
145     P(i,:) = P(i,:)/sumP(i);
146
147   end
148
149 end
150
151
152
153 % Check that the rowsums of Qc{1},...,Qc{r} are equal to 1
154
155
156
157 for i = 1:r
158
159   sumQi = sum(Qc{i}');
160
161   if sum(sumQi == 1) ~= length(Qc{i})
162
163     warning(['Rowsums of Q{' Zstr(i) '} not equal to 1. Renormalizing.']);
164
165     for j = 1:length(Qc{i})
166
167       Qc{i}(j,:) = Qc{i}(j,:)/sumQi(j);
168
169     end
170
171   end
172
173 end
174
175
176
177
178
179 % Make the transition matrix Q for the joint MC (X_k,Z_k)
180
181
182
183 Q = zeros(n*r,n*r);
184
185 I = 0:r:(n-1)*r;
186
187 for z=1:r
188
189   QQ = kron(Qc{z},P);
190
191   Q(I+z,:) = QQ(I+z,:);
192
193 end
194
195
196
197
198
199 % Stationary distribution = ro of Q
200
201
202
203 ro = mc2stat(Q);
204
205
206
207 % Generate random numbers
208
209
210
211 e=rand(T,1);
212
213
214
215 % Start values
216
217
218
219 e0 = e(1);
220
221 if isempty(init.z0) & isempty(init.x0)
222
223   x0z0 = min(find( e0<=cumsum(ro) ));
224
225   x0 = floor((x0z0+1)/r);
226
227   z0 = mod(x0z0-1,r)+1;
228
229 elseif isempty(init.x0)
230
231   z0 = init.z0;
232
233   rox0 = ro(z0:r:end); % Pick stat. distr. for regime z0
234
235   rox0 = rox0/sum(rox0);
236
237   x0 = min(find( e0<=cumsum(rox0) ));
238
239 elseif isempty(init.z0)
240
241   x0 = init.x0;
242
243   z0 = [];  % Start from stat. distr of P
244
245 else % Both z0 znd x0 are given
246
247   x0 = init.x0;
248
249   z0 = init.z0;
250
251 end
252
253
254
255 % Simulate Regime process
256
257
258
259 z = mcsim(P,T,z0);
260
261
262
263 % Simulate switching Markov chain
264
265
266
267 x=zeros(T,1);
268
269 x(1) = x0;   % First value
270
271
272
273 for k=2:T
274
275   Pi = Qc{z(k)};
276
277 %  eval(['Pi = P' Zstr(z(k)) ';']);
278
279   cumsumPi = cumsum(Pi')';
280
281   x(k) = min(find( e(k)<=cumsumPi(x(k-1),:) ));
282
283 end
284
285```

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