Changeset 4403

Show
Ignore:
Timestamp:
03/18/12 23:49:57 (6 years ago)
Author:
vdebuen
Message:

Refs #1136

Location:
tolp/OfficialTolArchiveNetwork/ArimaTools
Files:
6 added
4 modified

Legend:

Unmodified
Added
Removed
  • tolp/OfficialTolArchiveNetwork/ArimaTools/ArimaTools.tol

    r4397 r4403  
    3030#Embed "RandArima.tol"; 
    3131#Embed "arma_process.tol"; 
     32#Embed "OLS.tol"; 
     33#Embed "UnitRoot.tol"; 
    3234#Embed "Sampler.tol"; 
    3335#Embed "ARMA.Sampler.tol"; 
     
    5860 
    5961 
     62 
  • tolp/OfficialTolArchiveNetwork/ArimaTools/ARMA.Sampler.tol

    r4397 r4403  
    6767    Real _.L := Min(Max(5*_.s,_.p+_.q),_.N/4); 
    6868    Matrix _.acor := SubCol(AutoCor(difNoise, _.L),[[2]]); 
     69 
    6970    Set _.factor := For(1, _.numFactor, NameBlock(Real k) 
    7071    { 
     
    8081        Real L2 = Min(Max(5,p+q),L1,Floor(_.N/(4*s))); 
    8182        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) 
    8386      ]] 
    8487    }); 
     
    98101      Polyn ar = ChangeBDegree($factor.arima->AR,1/s); 
    99102      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) 
    101108    })) 
    102109  } 
     
    577584  }; 
    578585     
     586 
    579587  //////////////////////////////////////////////////////////////////////////// 
    580588  Set evalLogLikelihood_FastCholSea(Set arima) 
  • tolp/OfficialTolArchiveNetwork/ArimaTools/doc/Identificación baysesiana de modelos ARIMA multiestacionales.lyx

    r4399 r4403  
    127127\begin{array}{ccccc} 
    128128\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(1-B^{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} 
     129deg\left(\nabla\right)=D=\underset{i=1}{\overset{r}{\sum}}s_{r}d_{r} 
    130130\end{array} 
    131131\end{equation} 
     
    138138\begin{array}{ccccc} 
    139139\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} 
     140deg\left(\phi\right)=P=\underset{i=1}{\overset{r}{\sum}}s_{r}p_{r} 
    141141\end{array} 
    142142\end{equation} 
     
    150150\begin{array}{ccccc} 
    151151\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} 
     152deg\left(\theta\right)=Q=\underset{i=1}{\overset{r}{\sum}}s_{r}q_{r} 
    153153\end{array} 
    154154\end{equation} 
     
    250250 
    251251\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  
     252Aunque 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} 
     259d_{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} 
     270p_{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} 
     281q_{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 
     290El conjunto de todos los posibles vectores de grados ARIMA para una estacionalid 
     291ad 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 
     302y 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 
    258314\end_layout 
    259315 
    260316\begin_layout Standard 
    261317El 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} 
     328Pr\left[\left.\hat{\omega}\right|z\right]=\underset{\omega\in\Omega}{\max}\left\{ Pr\left[\left.\omega\right|z\right]\right\}  
     329\end{equation} 
     330 
     331\end_inset 
     332 
     333 
     334\end_layout 
     335 
     336\begin_layout Standard 
     337Si  
     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} 
     351Pr\left[\left.\omega\right|z\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 
     360La consideración misma del problema de determinar los grados ARIMA es tautológic 
     361amente 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} 
     372Pr\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 
     381Obviamente 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} 
     388d=\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 
     393y 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 
     405la serie de ruido diferenciado correspondiente es estacionaria y cuya probabilid 
     406ad condicionada por el ruido  
    279407\begin_inset Formula $z_{t}$ 
    280408\end_inset 
     
    284412 
    285413\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} 
     416Pr\left[\left.\omega\right|z\right]=Pr\left[\left.\varphi,d\right|z\right]=Pr\left[\left.\varphi\right|d,z\right]Pr\left[\left.\omega\right|z\right] 
     417\end{equation} 
     418 
     419\end_inset 
     420 
     421 
     422\end_layout 
     423 
     424\begin_layout Standard 
     425La consideración misma del problema de determinar los grados ARMA es tautológica 
     426mente equivalente a considerar diferentes alternativas de representación 
     427 como posibilidades prácticas. 
    299428 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} 
     438Pr\left[\varphi\left|d\right.\right]=Pr\left[\varphi'\left|d\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 
     447siempre y cuando el ruido diferenciado correspondiente a los grados de integraci 
     448ón 
    302449\begin_inset Formula  
    303450\begin{equation} 
  • tolp/OfficialTolArchiveNetwork/ArimaTools/Identify.tol

    r4397 r4403  
    55 
    66////////////////////////////////////////////////////////////////////////////// 
    7 Class @PolynDegreeSampler 
     7Class @PolDeg.Sampler 
    88////////////////////////////////////////////////////////////////////////////// 
    99{ 
     10  Static NameBlock Options =  
     11  [[ 
     12     NameBlock CandidateGen = [[ 
     13       Real MoveAll = 1; 
     14       Real MoveOne = 2; 
     15       Real MoveOneForced = 3 
     16     ]] 
     17  ]]; 
     18 
    1019  Real _.min; 
    1120  Real _.max; 
     
    1726 
    1827  //////////////////////////////////////////////////////////////////////////// 
    19   Real build_transition(Real lambda_)  
     28  Real build_transition(Real lambda_, Real candidateGen)  
    2029  //////////////////////////////////////////////////////////////////////////// 
    2130  {  
    2231    Real lambda := lambda_; 
     32     
    2333    Set _.transition := For(_.min, _.max, Set(Real x) 
    2434    { 
    2535      Set aux = For(_.min, _.max, Real(Real y) 
    2636      { 
    27       //If(x==y,0,lambda^Abs(x-y)) 
    28         lambda^Abs(x-y) 
     37        Case( 
     38          x!=y, lambda^Abs(x-y), 
     39          candidateGen==ArimaTools::@PolDeg.Sampler::Options:: 
     40          CandidateGen::MoveOneForced, 0, 
     41          1==1, 1) 
    2942      }); 
    3043      Real sum = SetSum(aux); 
     
    6073{ 
    6174  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); 
    7186    True  
    7287  } 
     
    93108  Set _.varMin = Copy(Empty); 
    94109  Set _.varMax = Copy(Empty); 
     110  Set _.varFree = Copy(Empty); 
    95111 
    96112  Set _.histogram = Copy(Empty); 
     
    107123    [[ 
    108124       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 
    110128    ]]; 
    111129    True 
     
    136154  //////////////////////////////////////////////////////////////////////////// 
    137155  { 
    138     Constant(_.varNum,1,0) 
     156    SetCol(_.varMin) 
    139157  }; 
    140158 
     
    153171      @FactorDegreeSampler fds = [[ 
    154172        Real period = dr->Periodicity; 
    155         @PolynDegreeSampler ar = [[ 
     173        @PolDeg.Sampler ar = [[ 
    156174          Real _.min = dr->MinARDegree, 
    157175          Real _.max = dr->MaxARDegree ]]; 
    158         @PolynDegreeSampler ma = [[ 
     176        @PolDeg.Sampler ma = [[ 
    159177          Real _.min = dr->MinMADegree, 
    160178          Real _.max = dr->MaxMADegree ]]; 
    161         @PolynDegreeSampler dif = [[ 
     179        @PolDeg.Sampler dif = [[ 
    162180          Real _.min = dr->MinDIFDegree, 
    163181          Real _.max = dr->MaxDIFDegree ]] 
     
    180198      SetOfReal(dr->MaxARDegree, dr->MaxMADegree, dr->MaxDIFDegree) 
    181199    })); 
     200    Set _.varFree := Select(Range(1,Card(_.varMin),1),Real(Real k) 
     201    { 
     202      _.varMin[k]<_.varMax[k] 
     203    }); 
    182204    Real _.varNum := Card(_.varName); 
    183205    Real stepSize := 0.5; 
     
    196218    Set EvalSet(_.degreeSampler,Real(@FactorDegreeSampler fds) 
    197219    { 
    198       fds::build_transition(stepSize) 
     220      fds::build_transition(stepSize,config::CandidateGen) 
    199221    }); 
    200222    True  
     
    257279      Real $ref::set.initialPoint($ref::get.initialPoint(?)); 
    258280/* * / 
     281      Set For(1,10,Real(Real iter) 
     282      { 
     283        Set arima = $ref::DrawStationary(?); 
     284        $ref::pi.logDens.sttCumAvr(arima) 
     285      }); 
     286/* * / 
    259287      Real $ref::optimizer.start(?); 
    260288      Set stopCriteria = @NameBlock(($($ref::_.optimizer):: 
    261289        methods[1])::stopCriteria); 
    262       Real $stopCriteria::maxTime := 1*(_.varNum^2); 
    263     //Real $stopCriteria::maxEval := 100; 
     290//    Real $stopCriteria::maxTime := 1*(_.varNum^2); 
     291      Real $stopCriteria::maxEval := 20; 
    264292      Real $stopCriteria::absoluteTolerance_target := .1; 
    265293      Real $stopCriteria::relativeTolerance_var := .01; 
     
    278306    SetCol(For(1,_.varNum,Real(Real k) 
    279307    { 
    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      }) 
    283314    })) 
    284315  }; 
     
    290321    SetSum(For(1,_.varNum,Real(Real k) 
    291322    { 
    292       @NameBlock dr = [[ _.degRef[k] ]]; 
    293323      Real m = MatDat(x,k,1); 
    294324      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      }) 
    296330    })) 
    297331  }; 
     
    303337  //////////////////////////////////////////////////////////////////////////// 
    304338  { 
    305     Real k = IntRand(1,_.varNum); 
     339    Real vf = IntRand(1,Card(_.varFree)); 
     340    Real k = _.varFree[vf]; 
    306341    Real _.lastMoving := k; 
    307342  //WriteLn("TRACE Q.draw.moveOne _.lastMoving = "<<_.lastMoving); 
     
    330365  //////////////////////////////////////////////////////////////////////////// 
    331366  { 
    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)) 
    333371  }; 
    334372 
     
    337375  //////////////////////////////////////////////////////////////////////////// 
    338376  { 
    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)) 
    340381  }; 
    341382 
     
    359400    Set arima = $arma::DrawStationary(?); 
    360401    $arma::pi.logDens.sttCumAvr(arima) 
    361   //$arma::pi.logDens.FactorACF(arima) 
    362   }; 
    363  
    364 /* */ 
     402  //$arma::pi.logDens($arma::arima2vector(arima)) 
     403  }; 
     404 
     405/* * / 
    365406  //////////////////////////////////////////////////////////////////////////// 
    366407  Real pi.logDens.last(Matrix x) 
     
    422463        Real i = MatDat(scsp,k,1); 
    423464        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) 
    427466        { 
    428467          Real deg = MatDat(_.mcmc,i,j);