Changeset 4403
 Timestamp:
 03/18/12 23:49:57 (6 years ago)
 Location:
 tolp/OfficialTolArchiveNetwork/ArimaTools
 Files:

 6 added
 4 modified
Legend:
 Unmodified
 Added
 Removed

tolp/OfficialTolArchiveNetwork/ArimaTools/ArimaTools.tol
r4397 r4403 30 30 #Embed "RandArima.tol"; 31 31 #Embed "arma_process.tol"; 32 #Embed "OLS.tol"; 33 #Embed "UnitRoot.tol"; 32 34 #Embed "Sampler.tol"; 33 35 #Embed "ARMA.Sampler.tol"; … … 58 60 59 61 62 
tolp/OfficialTolArchiveNetwork/ArimaTools/ARMA.Sampler.tol
r4397 r4403 67 67 Real _.L := Min(Max(5*_.s,_.p+_.q),_.N/4); 68 68 Matrix _.acor := SubCol(AutoCor(difNoise, _.L),[[2]]); 69 69 70 Set _.factor := For(1, _.numFactor, NameBlock(Real k) 70 71 { … … 80 81 Real L2 = Min(Max(5,p+q),L1,Floor(_.N/(4*s))); 81 82 Matrix acor = SubRow(_.acor,Range(s,Min(L2*s,_.N),s)); 82 Real L = Rows(acor) 83 Real L = Rows(acor); 84 Matrix cov = ARMA.ACF.Bartlett.Cov(acor,L,N); 85 Matrix cov_L = Choleski(cov) 83 86 ]] 84 87 }); … … 98 101 Polyn ar = ChangeBDegree($factor.arima>AR,1/s); 99 102 Polyn ma = ChangeBDegree($factor.arima>MA,1/s); 100 ARMA.ACF.Bartlett.LLH(ar,ma,$factor.acf::acor,$factor.acf::N) 103 Matrix tacov = ARMATACov(ar, ma, $factor.acf::L+1); 104 Matrix tacor = Sub(tacov,2,1,$factor.acf::L,1)*1/MatDat(tacov,1,1); 105 Matrix e = LTSolve(Tra($factor.acf::cov_L),tacor$factor.acf::acor); 106 0.5*Log(2*PI)0.5*MatSum(e$*e) 107 //ARMA.ACF.Bartlett.LLH(ar,ma,$factor.acf::acor,$factor.acf::N) 101 108 })) 102 109 } … … 577 584 }; 578 585 586 579 587 //////////////////////////////////////////////////////////////////////////// 580 588 Set evalLogLikelihood_FastCholSea(Set arima) 
tolp/OfficialTolArchiveNetwork/ArimaTools/doc/Identificación baysesiana de modelos ARIMA multiestacionales.lyx
r4399 r4403 127 127 \begin{array}{ccccc} 128 128 \nabla\left(B\right)=\underset{r=1}{\overset{m}{\prod}}\nabla_{r}\left(B^{s_{r}}\right) & \wedge & \nabla_{r}\left(B^{s_{r}}\right)=\left(1B^{s_{r}}\right)^{d_{r}} & \wedge & d_{r}\in\mathbb{N}\\ 129 deg\left(\nabla\right)= d=\underset{i=1}{\overset{r}{\sum}}s_{r}d_{r}129 deg\left(\nabla\right)=D=\underset{i=1}{\overset{r}{\sum}}s_{r}d_{r} 130 130 \end{array} 131 131 \end{equation} … … 138 138 \begin{array}{ccccc} 139 139 \phi\left(B\right)=\underset{r=1}{\overset{m}{\prod}}\phi_{r}\left(B^{s_{r}}\right) & \wedge & \phi_{r}\left(B^{s_{r}}\right)=1\underset{i=1}{\overset{p_{r}}{\sum}}\phi_{r,i}B^{s_{r}} & \wedge & p_{r}\in\mathbb{N}\\ 140 deg\left(\phi\right)= p=\underset{i=1}{\overset{r}{\sum}}s_{r}p_{r}140 deg\left(\phi\right)=P=\underset{i=1}{\overset{r}{\sum}}s_{r}p_{r} 141 141 \end{array} 142 142 \end{equation} … … 150 150 \begin{array}{ccccc} 151 151 \theta\left(B\right)=\underset{r=1}{\overset{m}{\prod}}\theta_{r}\left(B^{s_{r}}\right) & \wedge & \theta_{r}\left(B^{s_{r}}\right)=1\underset{i=1}{\overset{q_{r}}{\sum}}\theta_{r,i}B^{s_{r}} & \wedge & q_{r}\in\mathbb{N}\\ 152 deg\left(\theta\right)= q=\underset{i=1}{\overset{r}{\sum}}s_{r}q_{r}152 deg\left(\theta\right)=Q=\underset{i=1}{\overset{r}{\sum}}s_{r}q_{r} 153 153 \end{array} 154 154 \end{equation} … … 250 250 251 251 \begin_layout Standard 252 Aunque no es estrictamente necesario, supondremos que existe un grado mÃ¡ximo 253 254 \end_layout 255 256 \begin_layout Standard 257 Al conjunto de todos los posibles vectores de 252 Aunque no es estrictamente necesario para la argumentaciÃ³n, supondremos 253 que existe un grado mÃnimo y un grado mÃ¡ximo para cada componente 254 \end_layout 255 256 \begin_layout Standard 257 \begin_inset Formula 258 \begin{equation} 259 d_{r}\in\Omega_{r}^{\nabla}=\left\{ d_{r}^{},\cdots,d_{r}^{+}\right\} 260 \end{equation} 261 262 \end_inset 263 264 265 \end_layout 266 267 \begin_layout Standard 268 \begin_inset Formula 269 \begin{equation} 270 p_{r}\in\Omega_{r}^{\phi}=\left\{ p_{r}^{},\cdots,p_{r}^{+}\right\} 271 \end{equation} 272 273 \end_inset 274 275 276 \end_layout 277 278 \begin_layout Standard 279 \begin_inset Formula 280 \begin{equation} 281 q_{r}\in\Omega_{r}^{\theta}=\left\{ q_{r}^{},\cdots,q_{r}^{+}\right\} 282 \end{equation} 283 284 \end_inset 285 286 287 \end_layout 288 289 \begin_layout Standard 290 El conjunto de todos los posibles vectores de grados ARIMA para una estacionalid 291 ad concreta serÃ¡ 292 \end_layout 293 294 \begin_layout Standard 295 \begin_inset Formula 296 \begin{equation} 297 \Omega_{r}=\Omega_{r}^{\nabla}\times\Omega_{r}^{\phi}\times\Omega_{r}^{\theta}\subset\mathbb{N}^{3} 298 \end{equation} 299 300 \end_inset 301 302 y para el problema conjunto 303 \end_layout 304 305 \begin_layout Standard 306 \begin_inset Formula 307 \begin{equation} 308 \Omega=\underset{r=1}{\overset{m}{\prod}}\Omega_{r}\subset\mathbb{N}^{3m} 309 \end{equation} 310 311 \end_inset 312 313 258 314 \end_layout 259 315 260 316 \begin_layout Standard 261 317 El problema de la identificaciÃ³n un modelo ARIMA multiestacional consiste 262 en hallar los grados de los polinomios unitarios, autoregresivos y de media 263 mÃ³vil de cada una de las periodicidades 264 \end_layout 265 266 \begin_layout Standard 267 \begin_inset Formula 268 \begin{equation} 269 \Omega=\left\{ \Omega_{r}=\left\{ d_{r},p_{r},q_{r}\right\} \right\} _{r=1\ldots m} 270 \end{equation} 271 272 \end_inset 273 274 tales que la serie de ruido diferenciado 275 \begin_inset Formula $w_{t}$ 276 \end_inset 277 278 es estacionaria y cuya probabilidad condicionada por el ruido 318 en hallar el vector de enteros 319 \begin_inset Formula $\hat{\omega}\in\Omega$ 320 \end_inset 321 322 tales que su probabilidad condicionada por los datos sea mÃ¡xima 323 \end_layout 324 325 \begin_layout Standard 326 \begin_inset Formula 327 \begin{equation} 328 Pr\left[\left.\hat{\omega}\rightz\right]=\underset{\omega\in\Omega}{\max}\left\{ Pr\left[\left.\omega\rightz\right]\right\} 329 \end{equation} 330 331 \end_inset 332 333 334 \end_layout 335 336 \begin_layout Standard 337 Si 338 \begin_inset Formula $Pr\left[\omega\right]$ 339 \end_inset 340 341 es la probabilidad a priori de cada 342 \begin_inset Formula $\omega\in\Omega$ 343 \end_inset 344 345 entonces 346 \end_layout 347 348 \begin_layout Standard 349 \begin_inset Formula 350 \begin{equation} 351 Pr\left[\left.\omega\rightz\right]=Pr\left[\left.z\right\omega\right]Pr\left[\omega\right] 352 \end{equation} 353 354 \end_inset 355 356 357 \end_layout 358 359 \begin_layout Standard 360 La consideraciÃ³n misma del problema de determinar los grados ARIMA es tautolÃ³gic 361 amente equivalente a considerar diferentes alternativas de representaciÃ³n 362 como posibilidades prÃ¡cticas. 363 Por tanto, se tiene necesariamente que la probabilidad de cada combinaciÃ³n 364 de Ã³rdenes no puede ser nula a priori. 365 Si no hay mÃ¡s informaciÃ³n podemos suponer simplemente que todas las posibles 366 estructuras ARIMA son equiprobables a priori 367 \end_layout 368 369 \begin_layout Standard 370 \begin_inset Formula 371 \begin{equation} 372 Pr\left[\omega\right]=Pr\left[\omega'\right]\;\forall\omega,\omega'\in\Omega 373 \end{equation} 374 375 \end_inset 376 377 378 \end_layout 379 380 \begin_layout Standard 381 Obviamente el problema de la identificaciÃ³n ARIMA es separable en dos problemas 382 de identificaciÃ³n: identificar por un lado los Ã³rdenes de integraciÃ³n 383 \end_layout 384 385 \begin_layout Standard 386 \begin_inset Formula 387 \begin{equation} 388 d=\left(d_{1},\ldots,d_{m}\right)\in\Omega^{\nabla}=\underset{r=1}{\overset{m}{\prod}}\Omega_{r}^{\nabla}\subset\mathbb{N}^{m} 389 \end{equation} 390 391 \end_inset 392 393 y por otro los Ã³rdenes ARMA 394 \begin_inset Formula 395 \begin{equation} 396 \varphi=\left(p_{1},q_{1}\ldots,p_{m},q_{m}\right)\in\Omega^{\phi,\theta}=\underset{r=1}{\overset{m}{\prod}}\Omega_{r}^{\phi}\times\Omega_{r}^{\theta}\subset\mathbb{N}^{2m} 397 \end{equation} 398 399 \end_inset 400 401 402 \end_layout 403 404 \begin_layout Standard 405 la serie de ruido diferenciado correspondiente es estacionaria y cuya probabilid 406 ad condicionada por el ruido 279 407 \begin_inset Formula $z_{t}$ 280 408 \end_inset … … 284 412 285 413 \begin_layout Standard 286 La consideraciÃ³n misma del problema de determinar los grados ARMA 287 \end_layout 288 289 \begin_layout Standard 290 \begin_inset Formula 291 \begin{equation} 292 \Omega^{\phi,\theta}=\left\{ \Omega_{r}^{\phi,\theta}=\left\{ p_{r},q_{r}\right\} \right\} _{r=1\ldots m} 293 \end{equation} 294 295 \end_inset 296 297 es tautolÃ³gicamente equivalente a considerar diferentes alternativas de 298 representaciÃ³n como posibilidades prÃ¡cticas. 414 \begin_inset Formula 415 \begin{equation} 416 Pr\left[\left.\omega\rightz\right]=Pr\left[\left.\varphi,d\rightz\right]=Pr\left[\left.\varphi\rightd,z\right]Pr\left[\left.\omega\rightz\right] 417 \end{equation} 418 419 \end_inset 420 421 422 \end_layout 423 424 \begin_layout Standard 425 La consideraciÃ³n misma del problema de determinar los grados ARMA es tautolÃ³gica 426 mente equivalente a considerar diferentes alternativas de representaciÃ³n 427 como posibilidades prÃ¡cticas. 299 428 Por tanto, se tiene necesariamente que la probabilidad de cada combinaciÃ³n 300 de Ã³rdenes no puede ser nula, siempre y cuando el ruido diferenciado correspond 301 iente a los grados de integraciÃ³n 429 de Ã³rdenes no puede ser nula a priori. 430 Si no hay mÃ¡s informaciÃ³n podemos suponer simplemente que todas las posibles 431 estructuras ARMA condicionadas a una cierta estrura de integraciÃ³n son 432 equiprobables a priori 433 \end_layout 434 435 \begin_layout Standard 436 \begin_inset Formula 437 \begin{equation} 438 Pr\left[\varphi\leftd\right.\right]=Pr\left[\varphi'\leftd\right.\right]\;\forall\varphi,\varphi'\in\Omega^{\phi,\theta} 439 \end{equation} 440 441 \end_inset 442 443 444 \end_layout 445 446 \begin_layout Standard 447 siempre y cuando el ruido diferenciado correspondiente a los grados de integraci 448 Ã³n 302 449 \begin_inset Formula 303 450 \begin{equation} 
tolp/OfficialTolArchiveNetwork/ArimaTools/Identify.tol
r4397 r4403 5 5 6 6 ////////////////////////////////////////////////////////////////////////////// 7 Class @Pol ynDegreeSampler7 Class @PolDeg.Sampler 8 8 ////////////////////////////////////////////////////////////////////////////// 9 9 { 10 Static NameBlock Options = 11 [[ 12 NameBlock CandidateGen = [[ 13 Real MoveAll = 1; 14 Real MoveOne = 2; 15 Real MoveOneForced = 3 16 ]] 17 ]]; 18 10 19 Real _.min; 11 20 Real _.max; … … 17 26 18 27 //////////////////////////////////////////////////////////////////////////// 19 Real build_transition(Real lambda_ )28 Real build_transition(Real lambda_, Real candidateGen) 20 29 //////////////////////////////////////////////////////////////////////////// 21 30 { 22 31 Real lambda := lambda_; 32 23 33 Set _.transition := For(_.min, _.max, Set(Real x) 24 34 { 25 35 Set aux = For(_.min, _.max, Real(Real y) 26 36 { 27 //If(x==y,0,lambda^Abs(xy)) 28 lambda^Abs(xy) 37 Case( 38 x!=y, lambda^Abs(xy), 39 candidateGen==ArimaTools::@PolDeg.Sampler::Options:: 40 CandidateGen::MoveOneForced, 0, 41 1==1, 1) 29 42 }); 30 43 Real sum = SetSum(aux); … … 60 73 { 61 74 Real period; 62 @PolynDegreeSampler ar; 63 @PolynDegreeSampler ma; 64 @PolynDegreeSampler dif; 65 66 Real build_transition(Real lambda_) 67 { 68 Real ar::build_transition(lambda_); 69 Real ma::build_transition(lambda_); 70 Real dif::build_transition(lambda_); 75 @PolDeg.Sampler ar; 76 @PolDeg.Sampler ma; 77 @PolDeg.Sampler dif; 78 79 //////////////////////////////////////////////////////////////////////////// 80 Real build_transition(Real lambda_, Real candidateGen) 81 //////////////////////////////////////////////////////////////////////////// 82 { 83 Real ar::build_transition(lambda_, candidateGen); 84 Real ma::build_transition(lambda_, candidateGen); 85 Real dif::build_transition(lambda_, candidateGen); 71 86 True 72 87 } … … 93 108 Set _.varMin = Copy(Empty); 94 109 Set _.varMax = Copy(Empty); 110 Set _.varFree = Copy(Empty); 95 111 96 112 Set _.histogram = Copy(Empty); … … 107 123 [[ 108 124 Real Evaluator = 109 ArimaTools::@ARMA.Sampler::Options::Evaluator::Levinson 125 ArimaTools::@ARMA.Sampler::Options::Evaluator::Levinson; 126 Real CandidateGen = 127 ArimaTools::@PolDeg.Sampler::Options::CandidateGen::MoveAll 110 128 ]]; 111 129 True … … 136 154 //////////////////////////////////////////////////////////////////////////// 137 155 { 138 Constant(_.varNum,1,0)156 SetCol(_.varMin) 139 157 }; 140 158 … … 153 171 @FactorDegreeSampler fds = [[ 154 172 Real period = dr>Periodicity; 155 @Pol ynDegreeSampler ar = [[173 @PolDeg.Sampler ar = [[ 156 174 Real _.min = dr>MinARDegree, 157 175 Real _.max = dr>MaxARDegree ]]; 158 @Pol ynDegreeSampler ma = [[176 @PolDeg.Sampler ma = [[ 159 177 Real _.min = dr>MinMADegree, 160 178 Real _.max = dr>MaxMADegree ]]; 161 @Pol ynDegreeSampler dif = [[179 @PolDeg.Sampler dif = [[ 162 180 Real _.min = dr>MinDIFDegree, 163 181 Real _.max = dr>MaxDIFDegree ]] … … 180 198 SetOfReal(dr>MaxARDegree, dr>MaxMADegree, dr>MaxDIFDegree) 181 199 })); 200 Set _.varFree := Select(Range(1,Card(_.varMin),1),Real(Real k) 201 { 202 _.varMin[k]<_.varMax[k] 203 }); 182 204 Real _.varNum := Card(_.varName); 183 205 Real stepSize := 0.5; … … 196 218 Set EvalSet(_.degreeSampler,Real(@FactorDegreeSampler fds) 197 219 { 198 fds::build_transition(stepSize )220 fds::build_transition(stepSize,config::CandidateGen) 199 221 }); 200 222 True … … 257 279 Real $ref::set.initialPoint($ref::get.initialPoint(?)); 258 280 /* * / 281 Set For(1,10,Real(Real iter) 282 { 283 Set arima = $ref::DrawStationary(?); 284 $ref::pi.logDens.sttCumAvr(arima) 285 }); 286 /* * / 259 287 Real $ref::optimizer.start(?); 260 288 Set stopCriteria = @NameBlock(($($ref::_.optimizer):: 261 289 methods[1])::stopCriteria); 262 263 //Real $stopCriteria::maxEval := 100;290 // Real $stopCriteria::maxTime := 1*(_.varNum^2); 291 Real $stopCriteria::maxEval := 20; 264 292 Real $stopCriteria::absoluteTolerance_target := .1; 265 293 Real $stopCriteria::relativeTolerance_var := .01; … … 278 306 SetCol(For(1,_.varNum,Real(Real k) 279 307 { 280 @NameBlock dr = [[ _.degRef[k] ]]; 281 Real m = MatDat(x,k,1); 282 $dr::Q.draw(m) 308 If(_.varMin[k]==_.varMax[k], _.varMin[k], 309 { 310 @NameBlock dr = [[ _.degRef[k] ]]; 311 Real m = MatDat(x,k,1); 312 $dr::Q.draw(m) 313 }) 283 314 })) 284 315 }; … … 290 321 SetSum(For(1,_.varNum,Real(Real k) 291 322 { 292 @NameBlock dr = [[ _.degRef[k] ]];293 323 Real m = MatDat(x,k,1); 294 324 Real n = MatDat(y,k,1); 295 $dr::Q.logDens(m,n) 325 If(_.varMin[k]==_.varMax[k], If(m==n, 0, 1/0), 326 { 327 @NameBlock dr = [[ _.degRef[k] ]]; 328 $dr::Q.logDens(m,n) 329 }) 296 330 })) 297 331 }; … … 303 337 //////////////////////////////////////////////////////////////////////////// 304 338 { 305 Real k = IntRand(1,_.varNum); 339 Real vf = IntRand(1,Card(_.varFree)); 340 Real k = _.varFree[vf]; 306 341 Real _.lastMoving := k; 307 342 //WriteLn("TRACE Q.draw.moveOne _.lastMoving = "<<_.lastMoving); … … 330 365 //////////////////////////////////////////////////////////////////////////// 331 366 { 332 Q.draw.moveOne(x) 367 If(config::CandidateGen== 368 ArimaTools::@PolDeg.Sampler::Options::CandidateGen::MoveAll, 369 Q.draw.moveAll(x), 370 Q.draw.moveOne(x)) 333 371 }; 334 372 … … 337 375 //////////////////////////////////////////////////////////////////////////// 338 376 { 339 Q.logDens.moveOne(x,y) 377 If(config::CandidateGen== 378 ArimaTools::@PolDeg.Sampler::Options::CandidateGen::MoveAll, 379 Q.logDens.moveAll(x,y), 380 Q.logDens.moveOne(x,y)) 340 381 }; 341 382 … … 359 400 Set arima = $arma::DrawStationary(?); 360 401 $arma::pi.logDens.sttCumAvr(arima) 361 //$arma::pi.logDens .FactorACF(arima)362 }; 363 364 /* * /402 //$arma::pi.logDens($arma::arima2vector(arima)) 403 }; 404 405 /* * / 365 406 //////////////////////////////////////////////////////////////////////////// 366 407 Real pi.logDens.last(Matrix x) … … 422 463 Real i = MatDat(scsp,k,1); 423 464 Real frq = MatDat(scsp,k,2)/sampleLength; 424 Real cum := cum + frq; 425 Real cumFrq = Copy(cum); 426 [[frq,cumFrq]]<<For(1,_.varNum,Real(Real j) 465 [[frq]]<<For(1,_.varNum,Real(Real j) 427 466 { 428 467 Real deg = MatDat(_.mcmc,i,j);