" 64-bitars Firefox plugins
Övervakning av systemet "

Ett drama med 'R'

[Hoover_vote_1932.jpg]
Foto: Tony Misfit

En Mong de många intressanta applikationer som vi hittar i våra förråd av Ubuntu / Debian, är R , en statistisk paket med professionell anda ( jag använder Google, Merck, Pfizer och Bank of America ), som är resurser för nästan allt ni kan föreställa er, från grunderna, till de mest komplexa. Dess användning är inte, säger vi "mycket" vänlig "än kunskap om grunderna i S språket, är tolkningen av de olika funktionerna och dess många alternativ på något sätt omedelbart. Och även om du är intresserad av statistik, är ett bra paket för att rita data av något slag , och de mest fantasifulla sätt du kan tänka, och sparas via skript, en png, jpeg, pdf osv .... Men låt oss försöka det ... ett fall av mycket påtaglig.

Ärendet

Naturligtvis, om du läser detta från Spanien, är en av de särskiljande dragen i den ekonomiska krisen här den höga arbetslösheten vi har åstadkommit på några månader, och stor oro över den framtida utvecklingen av dessa uppgifter, de flesta att data är alla ett drama ... eller snarare många dramer ... Uppenbarligen är detta en arbetslöshet som är en följd och det skulle vara lämpligt att analysera i termer av andra variabler. Men eftersom detta är en ursäkt att leka med R vi kommer att behandla det i sig som en stokastisk process, och försöka härleda framtida utvecklingen med den typiska verktyg för tidsserieanalys av denna typ.

[Tasa_paro_sa_eurostat.png]
Källa och förberedelser: Eurostat

I forskning vi har en hel uppsättning av verktyg att hantera tidsserier och som inte gör analyser och prognoser av dem med ARIMA . Även om dessa modeller verkar atheoretical ibland lättsinnig, är i verktygslådan någon Econometrica, däribland Banco de España

Serien

Först måste vi välja källan för information. Detta i händelse av arbetslöshet är en grannlaga uppgift, eftersom det är beroende on som "ange" i statistiken ... Den bästa lösningen jag kan komma på är det månadsserier av Eurostat , den EPA may do bättre, men är kvartalsvis och vi skulle också sedan. Ett urval av cirka 80 observationer kommer att vara tillräcklig, så vi tar det från januari 2003 till maj 2009 (nominellt, justerat), det senast offentliggjorda. Den hämtas i csv med vad som direkt kan läsas av R. Värdena i kolumn INDICATOR_VALUE .

För det första att det är inte grafiken, få en uppfattning om hur serien: 

  2003 # strejk månad maj 2009 (Mediedelning EA15)
 yStart <- 2003
 DDIR <- "/ home / user / projekt / R / stopp / data /"
 # ES
 csv ( paste ( dDir , "Unemployment-NS-M-ES.csv" , sep = "" ) ) ESP <- läs. Csv (pasta (DDIR, "Arbetslöshet-NS-M-ES.csv", september = ""))
 pES [ 'INDICATOR_VALUE' ] , start = yStart , frequency = 12 ) paroES <- TS (PES ['INDICATOR_VALUE'], start = yStart, frekvens = 12)
 # EA15
 csv ( paste ( dDir , "/Unemployment-NS-M-EA15.csv" , sep = "" ) ) paroEA15 <- läs. csv (pasta (DDIR, "/ Arbetslöshet-NS-M-EA15.csv", september = ""))
 paroEA15 [ 'INDICATOR_VALUE' ] / 15 , start = yStart , frequency = 12 ) paroEA15 <- TS (paroEA15 ['INDICATOR_VALUE'] / 15, start = yStart, frekvens = 12)

 paroEA15 , paroES ) minmax <- intervall (paroEA15, paroES)
 X11 ()
 paroES , main = "Parados" , type = "l" , ylab = "x 1000 personas" , plot (paroES, huvud = "arbetslösa", type = "l", YLAB = "1000 personer x" 
             Enero 2003 - Mayo 2009" , XLab = "Månad \ n januari 2003 - maj 2009" 
             minmax [ 1 ] , minmax [ 2 ] + 1000 ) , ylim = c (minmax [1], [minmax 2] + 1000) 
             , sub = "(Fuente:Eurostat)" , lwd = 2 ) col = "red", sub = "(Källa: Eurostat)", lwd = 2)
 paroEA15 , col= "green4" , type = "l" , lwd = 2 ) linjer paroEA15 (, col = "green4" type = "l", lwd = 2)
 paroES ) [ 1 ] * dim ( paroES ) [ 2 ] upprepa <- dim (paroES) [1] * dim (paroES) [2]
 ts ( rep ( 1000 , repe ) , start = yStart , frequency = 12 ) , type = "l" , linjer (ts (rep 1000 (upprepade) start = yStart, frekvens = 12), type = "l" 
       , lty= "dashed" ) col = "grå", lty = "streckade")
 ts ( rep ( 4000 , repe ) , start = yStart , frequency = 12 ) , type = "l" , linjer (ts (rep 4000 (upprepade) start = yStart, frekvens = 12), type = "l" 
       , lty= "dashed" ) col = "grå", lty = "streckade")
 , legend= c ( "ES" , "EA15 (media)" ) , col= c ( 'red' , 'green4' ) , legend ("topright", legend = c ("ES", "EA15 (medel)"), col = c ("röda", "green4) 
        , lty= 1 , lwd= 1 ) BG = "gray99", lty = 1, lwd = 1)
 ( ) dev. off ()
I diagrammet ingår genomsnittet för de 15 länderna i euroområdet. Men denna jämförelse är inte helt korrekt (eller bättre något rätt) skulle ha att använda arbetslösheten för att detta ska vara betydande. Jag bara testar att representera flera serier i samma diagram.

[Paro_serie.png]

Model Selection: Box-Jenkins metod

Här kommer heuristisk del av processen. Ungefärliga värden av p, d, q.

Detta är ett icke-stationära serier, så vi måste differentiera och observera correlograms, och från dem härleda värdena för p och q. Men eftersom vi har en dator:) följer vi en metod som några diskreta uppgifter kriterier: AIC (Akaike) och / eller BIC (Bayes) . Närmar sig värdena för p och q (och d eftersom vi inte kommer att skilja i förväg) att minimera BIC och AIC.

AIC, kommer genom den mycket Arima i AIC del av klassen, och BIC är i paketet stats4 . Medan BIC stats4 Arima verkar inte fungera, så vi måste nöja oss med AIC, eller program som vi behöver kriterier. Eftersom det sistnämnda eftersom vi har den återstående variansen i komponenten klassen sigma2 Arima. Detta kan även inbegripa kriteriet Hannan-Quinn , HQIC.

Å andra sidan skulle behöva överväga säsongsvariation, med värden av p, d, q alternativ. Konsekvensen av cykliska mönster som visas i correlograms

Efter svettas lite med loopar, sluta vi att en modell SARIMA (0,2,0) x (0,1,1) kan ge oss goda resultat. 

  # ARIMA - Valda detaljhandel BIC och AIC
 arima ( pES [ 'INDICATOR_VALUE' ] , arbetslöshet. Arima <- Arima (PES ['INDICATOR_VALUE'] 
               0 , 2 , 0 ) , seasonal = list ( order = c ( 0 , 1 , 1 ) , period = 12 ) ) För = c (0, 2, 0), säsongens = lista (För = c (0, 1, 1), period = 12))
 arbetslöshet. Arima

 Kontrollera #
 c ( 1 , 1 ) ) pari (mfrow = c (1, 1))
 X11 ()
 arima ) # (test de Ljung-Box ERRONEO) tsdiag (unemployment. Arima) # (Ljung-Box test fel)
 ( ) dev. off ()

 # Kriterier för information ...
 # Paro.arima $ AIC
 paro . arima $ coef ) k = längd (unemployment. Arima $ koefficienten)
 pES ) T = längd (PES)
 ( paro . arima $ sigma2 ) + ( 2 * k ) / T #Akaike AIC <- log (unemployment. Arima $ sigma2) + (2 * k) / T # Akaike
 ( paro . arima $ sigma2 ) + ( k / T ) * log ( T ) #Bayesiano BIC <- log (unemployment. Arima $ sigma2) + (K / T) * log (T) # Bayesianska
 ( paro . arima $ sigma2 ) + ( ( 2 * k ) / T ) * log ( log ( T ) ) #Hannan-Quinn HQIC <- log (unemployment. Arima $ sigma2) + ((2 * k) / T) * log (log (T)) # Hannan-Quinn
 , format ( AIC ) , "BIC:" , format ( BIC ) , "HQIC:" , format ( HQIC ) ) Klistra in ("AIC", format (AIC), BIC: "format (BIC)," HQIC: "format (HQIC))

Och vi kan utvärdera autokorrelation av resthalter i diagrammet:

[Paro_arimadiag.png]

Här är några saker att notera första att korrelationen av eftersläpning 0, som alltid är 1. Detta är inte en bugg, bara i andra program bort nollan dröjsmål. Och för det andra att det är en bugg, är p-värden av Ljung-Box test om avfall, som returnerar tsdiag , som får oss att förkasta nollhypotesen. Av den anledningen inkluderar jag det i diagrammet.

Prediction

Nu kommer den intressanta delen med prognoserna, säger fram till slutet av året. 

  # Prediction
 predict ( paro . arima , n . ahead = 7 ) #Junio-Diciembre de 2009 paroES. förutsäga <- förutsäga (unemployment. Arima, n. framåt = 7) # juni till december, 2009

 # Plot prognos
 c ( 2009 , 6 ) #Comienzo predicción Pred. startar <- c (2009, 6) # start prognos
 predict $ pred + 2 * paroES . predict $ se UCS = paroES. Tippa $ Pred + 2 * paroES. Tippa $ är
 predict $ pred - 2 * paroES . predict $ se LCS = paroES. Tippa $ Pred - 2 * paroES. Tippa $ är
 predict <- min ( LCS , paroES ) min. paroES. förutsäga <- min (LCS, paroES)
 predict <- max ( UCS , paroES ) max. paroES. förutsäga <- max (UCS, paroES)
 ts ( UCS , start = pred . start , frequency = 12 ) ts. UCS <- TS (UCS, start = Pred. startar, frekvens = 12)
 ts ( LCS , start = pred . start , frequency = 12 ) ts. LCS <- TS (LCS, start = Pred. startar, frekvens = 12)

 X11 ()
 ( paroES , ts ( paroES . predict $ pred , start = pred . start , frequency = 12 ) , ts. plot (paroES, ts (paroES. förutsäga $ Pred, start = Pred. startar, frekvens = 12)
         Predicción" , type = "l" , viktigaste = "Arbetslösa SV \ n Prediction", type = "l"
         , xlab = "Mensual" , YLAB = "x 1000 människor" XLab = "Månadsvis" 
	 min . paroES . predict , max . paroES . predict ) , lwd = c ( 2 , 3 ) , ylim = c (min. paroES. förutsäga, max. paroES. förutsäga), lwd = c (2, 3)
	 "red" , "blue" ) , sub = "Enero 2003 - Diciembre 2009" ) col = c ("röd", "blå"), sub = "januari 2003 - december 2009)

 ts . UCS , col= "dark red" , lty= "dotted" , lwd = 2 ) linjer (ts. UCS, col = "mörkt röd", lty = "prickad", lwd = 2)
 ts . LCS , col= "dark green" , lty= "dotted" , lwd = 2 ) linjer (ts. LCS, col = "mörkgröna", lty = "prickad", lwd = 2) 
 ts ( rep ( 5000 , repe ) , start = yStart , frequency = 12 ) , type = "l" , linjer (ts (rep 5000 (upprepade) start = yStart, frekvens = 12), type = "l" 
       , lty= "dashed" ) col = "grå", lty = "streckade")
 ( ) dev. off ()

 # Skriv ut förutsägelse
 paroES . predict $ pred , start = pred . start , frequency = 12 ) Prediction <- TS (paroES. tippa $ Pred, start = Pred. Start, frekvens = 12)
 Predic , ts . UCS , ts . LCS ) resultat <- rbind (predikanter, ts. UCS, ts. LCS)
 <- c ( "Jun" , "Jul" , "Ago" , "Sep" , "Oct" , "Nov" , "Dic" ) colname (resultat) <- c ("juni", "juli", "sedan", "September", "Oktober" November "December")
 c ( "Predicción:" ) ) print (c ("Förutsägelse:"))
 resultat

Detta resultat ger oss: 

  "Prediction"
            Juni juli augusti september oktober november december
 Prediction 4436,84 4482.636 4592.294 4689.091 4834.444 4981.932 5102.101
 ts.UCS 4505,45 4636.053 4849.009 5064.884 5343.270 5636.430 5913.906
 ts.LCS 4368,23 4329.219 4335.578 4313.299 4325.618 4327.434 4290.297

Och vars graf:

[Paro_predict.png]

Ca 5 miljoner vid årets slut ... Jag hoppas inte det.

Jag upprepar att målet är en ursäkt för att illustrera användningen av R, inte en formell försök att förutsäga.

Taggar: ARIMA , nationalekonomi , statistik , Linux , R , Tidsserier

Denna artikel publicerades den Fredag 10 juli 2009 på 3:29 och är arkiverat under Program , nationalekonomi , statistik , Linux . Du kan följa någon svaren till denna post via RSS 2.0 . Du kan lämna en kommentar , eller göra trackback från din egen webbplats. Kommentarerna modereras och kommer inte att visas omedelbart.

Lämna en kommentar

Ange koden som visas: