Serie storica simulata su R
Qualcuno saprebbe linkarmi del materiale da cui "imparare" i comandi per simulare una serie storica con R?
Nel dettaglio, avrei bisogno di simulare una cosa del genere:
${(X_t\ |\ F_{t-1}:text{Poisson}(lambda_t); AAtinZZ) , (lambda_t=omega+alphaX_{t-1}):}$
(non mi viene la graffa grande su più righe
)
In pratica la distribuzione dell'osservazione al tempo $t$ condizionata alle osservazioni passate è una Poisson con un parametro che dipende dall'osservazione precedente. Si tratta di un INARCH(1) che può essere riscritto come: $X_t=omega+alphaX_{t-1}+epsilon_t$ con $epsilon_t=X_t-lambda_t$
Io avrei bisogno di simulare un certo numero di osservazioni da un processo del genere utilizzando R.
Nel dettaglio, avrei bisogno di simulare una cosa del genere:
${(X_t\ |\ F_{t-1}:text{Poisson}(lambda_t); AAtinZZ) , (lambda_t=omega+alphaX_{t-1}):}$
(non mi viene la graffa grande su più righe

In pratica la distribuzione dell'osservazione al tempo $t$ condizionata alle osservazioni passate è una Poisson con un parametro che dipende dall'osservazione precedente. Si tratta di un INARCH(1) che può essere riscritto come: $X_t=omega+alphaX_{t-1}+epsilon_t$ con $epsilon_t=X_t-lambda_t$
Io avrei bisogno di simulare un certo numero di osservazioni da un processo del genere utilizzando R.
Risposte
Aggiornamento =D
Consultando alcuni testi/articoli ho trovato il seguente codice che permette di simulare un processo di Poisson non omogeneo:
Il metodo consiste nello scegliere un $lambda>=lambda_t$, dove $lambda_t$ è la funzione d'intensità, che nell'esempio è la funzione seno mentre nel mio caso dovrebbe essere $omega+alpha*X_{t-1}$. E ogni $X_t$ dovrebbe essere un numero casuale estratto da una Poisson (esponenziale?) di parametro $lambda_t$ ($1/lambda_t$?).
Purtroppo ho fatto un po' di prove, ma non funziona
Tipo questa:
Consultando alcuni testi/articoli ho trovato il seguente codice che permette di simulare un processo di Poisson non omogeneo:
t<-0 I<-0 T<-100 lambda<-1.5 r<-runif(1) while(t<T){ t=t-1/lambda*log(r) if(r<sin(t)/lambda) I<-c(I,t) }
Il metodo consiste nello scegliere un $lambda>=lambda_t$, dove $lambda_t$ è la funzione d'intensità, che nell'esempio è la funzione seno mentre nel mio caso dovrebbe essere $omega+alpha*X_{t-1}$. E ogni $X_t$ dovrebbe essere un numero casuale estratto da una Poisson (esponenziale?) di parametro $lambda_t$ ($1/lambda_t$?).
Purtroppo ho fatto un po' di prove, ma non funziona

Tipo questa:
inarch<-function(omega,alfa,T,lambda){ I<-0 t<-0 x<-NULL lambda<-NULL x[0]<-0 for(t in 1:T){ lambda[t]<-omega+alfa*x[t-1] x[t]<-rpois(1,lambda[t]) } r<-runif(1) while(t<T){ t=t-1/lambda*log(r) if(r<=lambda[t]/lambda) I<-c(I,t) } }
Ciao,
secondo me può essere utile che copi i risultati sbagliati oppure gli errori del debugger, che ottieni.
Provo a sposare la discussione in Informatica, vediamo se lì qualcuno che conosce R sa aiutarti.
io cmq ti consiglierei di cambiare nome all'array lambda, non so come gestiche il typing R, di solito non è permesso l'overwriting nei linguaggi imperativi.
che tipo di proprietà deve avere la distribuzione degli eventi?
EDIT: domanda idiota, lo hai scritto in altra maniera sopra
secondo me può essere utile che copi i risultati sbagliati oppure gli errori del debugger, che ottieni.
Provo a sposare la discussione in Informatica, vediamo se lì qualcuno che conosce R sa aiutarti.
lambda[t]/lambda
io cmq ti consiglierei di cambiare nome all'array lambda, non so come gestiche il typing R, di solito non è permesso l'overwriting nei linguaggi imperativi.
"Arado90":
E ogni $X_t$ dovrebbe essere un numero casuale estratto da una Poisson (esponenziale?) di parametro $lambda_t$ ($1/lambda_t$?).
che tipo di proprietà deve avere la distribuzione degli eventi?
EDIT: domanda idiota, lo hai scritto in altra maniera sopra

io cmq ti consiglierei di cambiare nome all'array lambda, non so come gestiche il typing R, di solito non è permesso l'overwriting nei linguaggi imperativi.
C'hai ragione, correggiamo intanto questo!
inarch<-function(omega,alfa,T,lambda){ I<-0 t<-0 x<-NULL lam<-NULL x[0]<-0 for(t in 1:T){ lam[t]<-omega+alfa*x[t-1] x[t]<-rpois(1,lam[t]) } r<-runif(1) while(t<T){ t=t-1/lambda*log(r) if(r<=lam[t]/lambda) I<-c(I,t) } } inarch(1,0.7,1000,1.5) >Errore in rpois(1, lam[t]) : argomenti non validi
che tipo di proprietà deve avere la distribuzione degli eventi?
No, stavo pensando che $X_t$ è distribuita come una Poisson, però estraendo un solo numero casuale alla volta posso usare una esponenziale.
Il problema è che passando a rexp al posto di rpois il risultato finale è:
>There were 50 or more warnings (use warnings() to see the first 50)
xD
Guarda io di processi stocastici non ne so molto, tuttavia scrivere:
genera un errore, mentre scrivere:
credo faccia quanto voluto. Semplicemente i vettori in R cominciano dall'elemento 1 e non dall'elemento 0. Basta posticipare il ciclo di 1 e dovrebbe funzionare.
> x<-NULL > x[0]<-0 > x numeric(0)
genera un errore, mentre scrivere:
x<-NULL > x[1]<-0 > x [1] 0
credo faccia quanto voluto. Semplicemente i vettori in R cominciano dall'elemento 1 e non dall'elemento 0. Basta posticipare il ciclo di 1 e dovrebbe funzionare.
Dannazione, grazie mille!! Mi era sfuggito questo particolare!
Adesso il programma quantomeno funziona =D
Resta da verificare la correttezza di quanto sto facendo, cioè se la serie di numeri che finalmente il programma restituisce è del tipo che serve a me.
Adesso il programma quantomeno funziona =D
Resta da verificare la correttezza di quanto sto facendo, cioè se la serie di numeri che finalmente il programma restituisce è del tipo che serve a me.