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 validiche 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.