Come regredire più serie su serie singole condizionate su variabili di raggruppamento?

Question

Ritengo che il mio problema di base sia come regredire più serie su una singola serie. Anche se le mie serie non sono uguali nel tempo, anche quando uso serie di durata pari per lo stock e il benchmark (posso fornire dati che ho reso uguali manualmente se necessario) ottengo un errore. Voglio stimare un modello di mercato (cioè, regredire i rendimenti di un titolo sul rendimento del benchmark per ogni giorno, per tutti gli stock) e fare un quadro di dati di valori beta dalla regressione in formato lungo. Quindi, per il campione fornito, ci saranno 4 valori beta (2 per ABC e 2 per XYZ) nel riquadro dati del valore beta. Questo è un esempio di due corsi azionari

idf <- structure(list(Firm = c("ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", 
    "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "ABC", "XYZ", "XYZ", "XYZ", 
    "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", "XYZ", 
    "XYZ"), Date = structure(c(NA, 1451642400, 1451646000, 1451649600, 1451653200, 
    1451656800, 1451660400, 1451664000, 1451898000, 1451901600, 1451905200, 
    1451908800, 1451912400, 1451916000, 1451919600, NA, 1451642400, 1451646000, 
    1451649600, 1451653200, 1451656800, 1451660400, 1451664000, 1451898000, 
    1451901600, 1451905200, 1451908800, 1451912400, 1451916000, 1451919600), 
    tzone = "UTC", class = c("POSIXct", "POSIXt")), Price = c(1270.9, 1277, 
    1273.25, 1273.85, 1273.75, 1272, 1265.35, 1265.35, 1253.1, 1248.1, 1242, 
    1248.15, 1241.1, 1246.5, 1242.5, 225.75, 225.7, 225.5, 225.45, 228.6, 227.7, 
    227.8, 227.8, 226, 225.1, 222.35, 222.25, 221.1, 221.2, 220.7), rt = c(NA, 
    0.00478826614451489, -0.0029408902678254, 0.000471124032113579, 
    -7.8505259892836e-05, -0.00137484063686699, -0.00524170116535849, 0, 
    -0.00972828256347036, -0.00399808624683118, -0.00489941143098349, 
    0.00493947152112462, -0.00566437187907365, 0.00434154082813709, 
    -0.00321414499282824, NA, -0.000221508473604359, -0.000886524880757023, 
    -0.000221754075639957, 0.0138753464936165, -0.00394477829101625, 
    0.000439077943387822, 0, -0.00793305174079517, -0.00399025135959974, 
    -0.0122920309559813, -0.000449842562691316, -0.00518778653061336, 
    0.000452181784779349, -0.00226295638548901), day = structure(c(NA, 16801, 
    16801, 16801, 16801, 16801, 16801, 16801, 16804, 16804, 16804, 16804, 16804, 
    16804, 16804, NA, 16801, 16801, 16801, 16801, 16801, 16801, 16801, 16804, 
    16804, 16804, 16804, 16804, 16804, 16804), class = "Date")), 
    .Names = c("Firm", "Date", "Price", "rt", "day"), 
    row.names = c(NA, -30L), class = c("tbl_df", "data.frame")) 

e un campione di indice di riferimento

imdf <- structure(list(Date = structure(c(NA, 1451642400, 1451646000, 1451649600, 
    1451653200, 1451656800, 1451660400, 1451664000, 1451898000, 1451901600, 
    1451905200, 1451908800, 1451912400, 1451916000, 1451919600, 1451923200), 
    class = c("POSIXct", "POSIXt"), tzone = "UTC"), Price = c(3443.1, 3450.5, 
    3453.85, 3447.9, 3456.9, 3468.45, 3472.1, 3472.1, 3484.75, 3466.45, 3417.5, 
    3416.05, 3401, 3425.75, 3425.25, 3425.25), rt = c(NA, 0.0021469197059254, 
    0.000970402793278424, -0.00172420081111291, 0.00260688364525663, 
    0.00333557458000655, 0.00105178994070698, 0, 0.0036367074012027, 
    -0.00526529010184795, -0.0142217259100423, -0.000424376794422088, 
    -0.00441540679648789, 0.00725091981911596, -0.000145964093093198, 
    0), day = structure(c(NA, 16801, 16801, 16801, 16801, 16801, 16801, 16801, 
    16804, 16804, 16804, 16804, 16804, 16804, 16804, 16804), class = "Date")), 
    .Names = c("Date", "Price", "rt", "day"), row.names = c(NA, -16L), 
    class = c("tbl_df", "tbl", "data.frame")) 

seguente è il mio codice di base dplyr per la stima attuale beta per gli stock ABC e XYZ per due giorni regredendo loro intraday ritorno sul rendimento del benchmark.

require(dplyr) 
q3 <- idf %>% 
    group_by(Firm, day) %>% 
    summarise(Beta = PerformanceAnalytics::CAPM.beta(idf$rt, imdf$rt)) %>% 
    na.omit() 

Ma dà il seguente errore:

Error: The data cannot be converted into a time series. If you are trying to pass in names from a data object with one column, you should use the form 'data[rows, columns, drop = FALSE]'. Rownames should have standard date formats, such as '1985-03-15'. 

Per

q3 <- idf %>% 
group_by(Firm, day) %>% 
summarise(Beta = coef(summary(lm(idf$rt~imdf$rt)))[2,1]) %>% 
na.omit() 

ottengo il seguente errore:

Error: error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = idf$rt ~ imdf$rt, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt')` 

Ho anche provato la seguente basata sul complesso funzione:

beta.est= function(x) coef(summary(lm(x~imdf$rt)))[2,1] 
betas = aggregate(cbind(rt) ~ Firm + day , idf, FUN = beta.est) 

Ma dà anche lo stesso errore:

Error in summary(lm(x ~ imdf$rt)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = x ~ imdf$rt, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt') 

Poi ho provato un 'ciclo':

betas=list() 
for(i in idf$Firm){ 
for(j in idf$day){ 
    betas$day= coef(summary(lm(idf[i,4]|j~imdf$rt|j)))[2,1]}} 

Ancora una volta, ottengo questo:

Error in summary(lm(idf[i, 4] | j ~ imdf$rt | j)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = idf[i, 4] | j ~ imdf$rt | j, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt | j') 

Quando uso serie di durata pari per lo stock e il benchmark, ottengo il seguente errore r:

Error in summary(lm(idf[i, 4] | j ~ imdf$rt | j)) : error in evaluating the argument 'object' in selecting a method for function 'summary': Error in model.frame.default(formula = idf[i, 4] | j ~ imdf$rt | j, drop.unused.levels = TRUE) : variable lengths differ (found for 'imdf$rt | j') 

Sebbene i codici di cui sopra sono per stimare valori di beta correnti, sarebbe bello se posso anche regredire ritardo di indice di riferimento in cantiere. Gentilmente aiutami.

Answer 1

Spero non ti dispiaccia che utilizzo alcuni pacchetti aggiuntivi.Ho il mio parere una combinazione di dplyr, tidyr, purrr e broom davvero simplifys lavorare con più regressioni lineari:

per assicurarci di non confondiamo i nomi delle colonne di dati con i nomi delle colonne di riferimento abbiamo semplicemente aggiungere il nome di riferimento:

names(imdf) <- paste("bm1", names(imdf), sep ="_") 

È inoltre possibile aggiungere una variabile Data lag e/o partecipare a un altro benchmark

set.seed(123) 
imdf$bm1_Date_lag <- imdf$bm1_Date - 3600  # shift Date by 1 hour 
imdf$bm2_rt <- rnorm(nrow(imdf), sd = 0.005) # add a dummy benchmark 

Prima di tutto bisogna raggruppare i dati per Firm e day, unitevi ai dati di riferimento su Date e annidatelo, questo crea per ogni gruppo un data.frame chiamato data contenente tutte le altre variabili.

idf_grouped <- idf %>% group_by(Firm, day) %>% 
    left_join(imdf, by = c("Date" = "bm1_Date")) %>% 
    nest() 

Rimuovere i gruppi senza day valori in modo che siano in grado di adattarsi a un modello lineare su ciascun gruppo e aggiungerlo come nuova colonna con map da purrr. Puoi farlo per più modelli in un solo passaggio.

idf_grouped <- idf_grouped %>% filter(!is.na(day)) %>% 
    mutate(model = map(data,~lm(rt~bm1_rt, data = .)), 
      model2 = map(data,~lm(rt~bm2_rt, data = .))) 

per estrarre i coefficienti di modelli lineari che preferisco tidy da broom perché restituisce un data.frame che è - nel contesto di dplyr - più facile da gestire. unnest() trasforma i dati in un unico data.frame. Perché non si è interessati al intercetta filtriamo per bm1_rt

idf_grouped %>% 
    unnest(bm1 = model %>% map(tidy, quick = TRUE), 
      bm2 = model2 %>% map(tidy, quick = TRUE),.sep = "_") %>% 
    filter(bm1_term == "bm1_rt") 
## Firm  day bm1_term bm1_estimate bm2_term bm2_estimate 
## <chr>  <date> <fctr>  <dbl> <fctr>  <dbl> 
## 1 ABC 2016-01-01 bm1_rt 0.08879514 bm2_rt -0.2781275 
## 2 ABC 2016-01-04 bm1_rt 0.25888489 bm2_rt 0.3845765 
## 3 XYZ 2016-01-01 bm1_rt 0.66791986 bm2_rt -0.2891714 
## 4 XYZ 2016-01-04 bm1_rt 0.45735812 bm2_rt -0.5014824 

Per calcolare questi beta su lag benchmark appena si uniscono per "Date" = "bm1_Date_lag" durante la prima fase.

Modifica: Per implementare un reso di settore è necessario disporre di una mappatura a cui appartiene ciascuna azienda. Per scopo di esempio ho semplicemente aggiunto un'altra impresa "DEF" come una copia di "XYZ" che traccio per lo stesso settore come "ABC"

idf_2 <- idf %>% filter(Firm == "XYZ") %>% mutate(Firm = "DEF") %>% bind_rows(idf) 
firm_map <- data.frame(Firm = c("ABC", "DEF", "XYZ"), Industry = c(1,1,2), stringsAsFactors = FALSE) 

semplicemente aderire a questo in idf_2

idf_map <- idf_2 %>% left_join(firm_map, by = "Firm") 

e calcolare ad esempio. rendimento medio di ogni settore

idf_map %>% left_join(idf_map %>% group_by(Industry, Date) %>% 
    summarise(ind_rt = mean(rt, na.rm = TRUE)), by = c("Industry", "Date")) 

Ora ind_rt può essere usato come variabile esplicativa nella regressione.