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: This function is called by:

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