Stochastic portfolio optimisation and out-of-sample testing using R

Author: Andreï Victorovitch Kostyrka, University of Luxembourg.

Last updated: 2018-07-31.

One might ask a question, “What is it like to try to pick an optimal portfolio and invest in stocks?” Well, here I have prepared some R code that does the following thing:

  • Downloads stock data from Yahoo Finance from 2014 for a selection of stocks;
  • Puts all stocks together across the maximum set of available dates and nicely imputes missing values using Kalman filter with structural modelling;
  • Tries all possible combinations or portfolios (as many as possible, using Halton pseudo-random low-discrepancy sequences, doing many iterations and saving the best result from every iteration to save memory);
  • Creates amazing colourful plots for subsequent analysis;
  • Returns statistics for equally weighted portfolios, max-return under equally weighted risk, min-risk under equally weighted return, and maximum Sharpe ratio (for simplicity, the ratio of returns to risk, assuming there is no risk-free asset);
  • Tests these 4 portfolios out of sample for 1 year (“if one started trading one year ago, what would one have obtained”).

Two sets of stocks were chosen: most lucrative U.S. companies and most lucrative Russian companies. This is what one gets if one chooses 13 top U.S. companies:

library(BatchGetSymbols)
library(zoo)
library(imputeTS) # For imputation
library(randtoolbox) # For Sobel sampling
library(parallel)
 
# http://fortune.com/2017/06/07/fortune-500-companies-profit-apple-berkshire-hathaway/
# https://www.swiftpassportservices.com/blog/profitable-companies-russia/
 
syms <- c("GOOG", "AAPL", "JPM", "BRK-A", "WFC", "MCD", "BAC", "MSFT", "AMZN", "JNJ", "C", "MO", "GE")
symnames <- c("Google", "Apple", "J.P. Morgan", "Berkshire Hathaway", "Wells Fargo", "McDonald's",
              "Bank of America", "Microsoft", "Amazon", "Johnson and Johnson", "Citigroup", "Altria group", "General Electric")
cat("Calculating optimal portfolio for US stocks\n")
 
# syms <- c("SGTZY", "ROSN.ME", "OGZPY", "LUKOY", "SBRCY", "TRNFP.ME", "NILSY", "ALRS.ME", "YNDX")
# symnames <- c("SurgutNefteGas", "RosNeft", "GazProm", "LukOil", "SberBank", "TransNeft",
#               "Norilsk Nickel", "Alrosa", "Yandex")
# cat("Calculating optimal portfolio for Russian stocks\n")
 
first.date <- as.Date("2014-01-01")
last.date <- as.Date("2018-08-01")
end.ins <- as.Date("2017-08-01")
 
dat <- BatchGetSymbols(tickers = syms, first.date = first.date, last.date = last.date)
 
allt <- dat$df.tickers
allt <- allt[!is.na(allt$price.adjusted), ]
dates <- sort(unique(allt$ref.date))
prices <- data.frame(dates)
prices[, 2:(length(syms)+1)] <- NA
colnames(prices)[2:(length(syms)+1)] <- syms
for (i in syms) {
  datesi <- allt$ref.date[as.character(allt$ticker)==i]
  pricesi <- allt$price.adjusted[as.character(allt$ticker)==i]
  prices[match(datesi, dates), i] <- pricesi
  prices[, i] <- zoo(prices[, i], order.by = prices$dates)
  prices[, i] <- na.kalman(prices[, i])
}
prices2 <- prices
prices[,1] <- NULL
 
co <- rainbow(length(syms), end=0.75)
png("1-relative-prices.png", 700, 400, type="cairo")
plot(prices[,1]/mean(prices[,1]), col=co[1], lwd=1.5, main="Prices divided by their in-sample mean", xlab="Time", ylab="Ratio")
for (i in 2:length(syms)) {
  lines(prices[,i]/mean(prices[,i]), col=co[i])
}
legend("topleft", syms, col=co, lwd=1)
dev.off()
 
rs <- data.frame(rep(NA, nrow(prices)-1))
for (i in syms) {
  rs[, i] <- diff(log(prices[, i]))
}
rs[,1] <- NULL
 
w0 <- rep(1/length(syms), length(syms))
names(w0) <- syms
rmean <- zoo(as.numeric(as.matrix(rs) %*% w0), order.by = time(rs[,1]))
 
a <- table(as.numeric(format(time(rmean), "%Y")))
daysperyear <- mean(a[1:(length(a)-1)]) # Because the last year is incomplete
 
rsfull <- rs # Full sample
rs <- rs[time(rs[,1]) <= end.ins, ] # Training sample
rsnew <- rsfull[time(rsfull[,1]) > end.ins, ] # Test sample
 
meanrs <- colMeans(rs)
covrs <- cov(rs)
 
musigma <- function(w) {
  mu <- sum(meanrs * w)*daysperyear
  s <- sqrt((t(w) %*% covrs %*% w) * daysperyear)
  return(c(mu=mu, sigma=s, ratio=mu/s))
}
 
m0 <- musigma(w0)
 
num.workers <- if (.Platform$OS.type=="windows") 1 else detectCores()
MC <- 108000
iters <- 100
cc <- 150 # How many colours
cols <- rainbow(cc, end = 0.75)
sharpecuts <- seq(-0.1, 1.5, length.out = cc)
 
png("2-portfolios.png", 700, 600, type="cairo")
 
# First iteration
cat("Calculating", MC, "portfolios for the first time.\n")
ws <- halton(MC, length(syms), init=TRUE)
ws <- ws/rowSums(ws)
colnames(ws) <- syms
res <- mclapply(1:MC, function(i) musigma(ws[i, ]), mc.cores = num.workers)
res <- as.data.frame(matrix(unlist(res, use.names = FALSE), nrow=MC, byrow = TRUE))
colnames(res) <- c("mean", "sd", "ratio")
prcmu0 <- ecdf(res$mean)(m0[1]) # Which quantile is the mean
prcsd0 <- ecdf(res$sd)(m0[2]) # Which quantile is the sd
qmplus <- quantile(res$mean, prcmu0+0.005)
qmminus <- quantile(res$mean, prcmu0-0.005)
qsplus <- quantile(res$sd, prcsd0+0.005)
qsminus <- quantile(res$sd, prcsd0-0.005)
meanvicinity <- which(res$mean < qmplus & res$mean > qmminus)
sdvicinity <- which(res$sd < qsplus & res$sd > qsminus)
 
plot(res[,2:1], pch=16, col=cols[as.numeric(cut(res$ratio, breaks = sharpecuts))],
     xlim=range(res$sd)+c(-0.03, 0.03), ylim=pmax(range(res$mean)+c(-0.03, 0.03), c(0, -10)),
     main="Risk and return for various portfolios", xlab="Annualised volatility", ylab="Annualised log-return")
 
sharpes <- mus <- sds <- array(NA, dim=c(iters+1, length(syms)))
sharpes[1, ] <- ws[which.max(res$ratio), ]
mus[1, ] <- ws[which(res$mean == max(res$mean[sdvicinity])), ]
sds[1, ] <- ws[which(res$sd == min(res$sd[meanvicinity])), ]
pairsm <- pairss <- array(NA, dim=c(iters+1, 3))
pairsm[1, ] <- res$mean[c(which.max(res$ratio), which(res$mean == max(res$mean[sdvicinity])), which(res$sd == min(res$sd[meanvicinity])))]
pairss[1, ] <- res$sd[c(which.max(res$ratio), which(res$mean == max(res$mean[sdvicinity])), which(res$sd == min(res$sd[meanvicinity])))]
 
for (i in 1:iters) {
  cat("Calculating", MC, "portfolios, iteration", i, "out of", iters, "\n")
  # ws <- sobol(MC, length(syms), scrambling = 1, init=FALSE) # There is a bug in randtoolbox: it starts returning numbers outside [0, 1]
  ws <- halton(MC, length(syms), init=FALSE)
  ws <- ws/rowSums(ws)
  colnames(ws) <- syms
  res <- mclapply(1:MC, function(i) musigma(ws[i, ]), mc.cores = num.workers)
  res <- as.data.frame(matrix(unlist(res, use.names = FALSE), nrow=MC, byrow = TRUE))
  colnames(res) <- c("mean", "sd", "ratio")
  meanvicinity <- which(res$mean < qmplus & res$mean > qmminus)
  sdvicinity <- which(res$sd < qsplus & res$sd > qsminus)
  maxsh <- which.max(res$ratio)
  maxm <- which(res$mean == max(res$mean[sdvicinity]))
  minsd <- which(res$sd == min(res$sd[meanvicinity]))
  sharpes[i+1, ] <- ws[maxsh, ]
  mus[i+1, ] <- ws[maxm, ]
  sds[i+1, ] <- ws[minsd, ]
  pairsm[i+1, ] <- res$mean[c(maxsh, maxm, minsd)]
  pairss[i+1, ] <- res$sd[c(maxsh, maxm, minsd)]
  points(res[,2:1], pch=16, col=cols[as.numeric(cut(res$ratio, breaks = sharpecuts))])
  rm(res); gc()
}
 
points(as.numeric(pairss), as.numeric(pairsm))
 
sharpems <- apply(sharpes, 1, musigma)
wmaxsharpe <- sharpes[which.max(sharpems[3,]), ]
meanms <- apply(mus, 1, musigma)
wmaxmean <- mus[which.max(meanms[1,]), ] # Maximum returns while SD is the same
sdms <- apply(sds, 1, musigma)
wminsd <- sds[which.min(sdms[2,]), ] # Minimum SD while mean is the same
 
m1 <- musigma(wmaxsharpe)
m2 <- musigma(wmaxmean)
m3 <- musigma(wminsd)
finalstats <- cbind(rbind(w0, wmaxsharpe, wmaxmean, wminsd), rbind(m0, m1, m2, m3))
rownames(finalstats) <- c("Equal", "MaxSharpe", "MaxReturn", "MinRisk")
colnames(finalstats)[(length(syms)+1):(length(syms)+3)] <- c("RInS", "VInS", "SRInS")
 
 
points(m0[2], m0[1], cex=3, pch=16)
points(c(m1[2], m2[2], m3[2]), c(m1[1], m2[1], m3[1]), cex=2, pch=16)
dev.off()
 
# And now checking these revenues on the test sample
r0 <- zoo(as.numeric(as.matrix(rsnew) %*% w0), order.by = time(rsnew[,1]))
r1 <- zoo(as.numeric(as.matrix(rsnew) %*% wmaxsharpe), order.by = time(rsnew[,1]))
r2 <- zoo(as.numeric(as.matrix(rsnew) %*% wmaxmean), order.by = time(rsnew[,1]))
r3 <- zoo(as.numeric(as.matrix(rsnew) %*% wminsd), order.by = time(rsnew[,1]))
 
png("3-test-sample-returns.png", 700, 400, type="cairo")
plot(r0, main="Out-of-sample returns of weighted portfolios", ylab="Return", xlab="Time")
lines(r1, col=2)
lines(r2, col=3)
lines(r3, col=4)
legend("topright", c("Equal", "Sharpe", "Max. return", "Min. volatility"), col=1:4, lwd=1)
dev.off()
 
png("4-test-sample.png", 700, 400, type="cairo")
plot(exp(cumsum(r0)), ylim=c(0.6, 1.4), main="Portfolio growth over the last year in the test sample", ylab="Relative growth", xlab="Time")
abline(h=1)
lines(exp(cumsum(r1)), col=2)
lines(exp(cumsum(r2)), col=3)
lines(exp(cumsum(r3)), col=4)
legend("bottomleft", c("Equal", "Sharpe", "Max. return", "Min. volatility"), col=1:4, lwd=1)
dev.off()
 
finals <- cbind(r0, r1, r2, r3)
meanf <- colMeans(finals)*daysperyear
sdf <- apply(finals, 2, sd)*sqrt(daysperyear)
finalsoos <- cbind(meanf, sdf, meanf/sdf)
colnames(finalsoos) <- c("ROOS", "VOOS", "SRROS")
finalstats <- cbind(finalstats, finalsoos)
colnames(finalstats)[1:length(syms)] <- symnames
 
library(stargazer)
stargazer(t(finalstats), summary = FALSE, type="html", out="a.html")
 
write.csv(prices2, "stocks.csv", row.names = FALSE)

This is how the 13 U.S. stocks behaved.

Relative prices of top U.S. stocks

This is how 10 million random portfolios behaved (black dots indicate 100 best portfolios under the three criteria; big black dot for the equally weighted portfolio).

Return and risk of various U.S. portfolios

These are portfolio weights for U.S. companies with their performance.

Equal MaxSharpe MaxReturn MinRisk
Google 0.077 0.011 0.004 0.052
Apple 0.077 0.132 0.128 0.105
J.P. Morgan 0.077 0.133 0.099 0.028
Berkshire Hathaway 0.077 0.019 0.040 0.170
Wells Fargo 0.077 0.003 0.008 0.025
McDonald's 0.077 0.164 0.148 0.192
Bank of America 0.077 0.030 0.034 0.036
Microsoft 0.077 0.078 0.180 0.022
Amazon 0.077 0.088 0.161 0.005
Johnson & Johnson 0.077 0.138 0.064 0.138
Citigroup 0.077 0.003 0.003 0.009
Altria group 0.077 0.198 0.131 0.142
General Electric 0.077 0.003 0.001 0.078
Return in-sample 0.145 0.179 0.190 0.145
Risk in-sample 0.142 0.128 0.141 0.119
Sharpe ratio in-sample 1.023 1.395 1.345 1.217
Return in 2017.07–2018.06 0.135 0.163 0.241 0.053
Risk in 2017.07–2018.06 0.144 0.134 0.146 0.130
Sharpe ratio in 2017.07–2018.06 0.941 1.218 1.657 0.404

This is how these four U.S. portfolios look out of sample.

Returns for different optimal U.S. portfolios out-of-sample

If one had invested in these four portfolios a year before, this is what one would have harvested by now.

Cumulative growth of one's investments into U.S. stocks


And this is the result for top 9 Russian companies.

This is how the 9 Russian stocks behaved.

Relative prices of top Russian stocks

This is how 10 million random portfolios behaved (black dots indicate 100 best portfolios under the three criteria; big black dot for the equally weighted portfolio).

Return and risk of various Russian portfolios

Note that some portfolios can actually yield negative returns, and the volatility is generally higher!

These are portfolio weights for Russian companies with their performance.

Equal MaxSharpe MaxReturn MinRisk
SurgutNefteGas 0.111 0.005 0.009 0.175
RosNeft 0.111 0.162 0.031 0.269
GazProm 0.111 0.012 0.068 0.067
LukOil 0.111 0.018 0.083 0.037
SberBank 0.111 0.038 0.031 0.014
TransNeft 0.111 0.248 0.016 0.213
Norilsk Nickel 0.111 0.079 0.381 0.044
Alrosa 0.111 0.413 0.358 0.049
Yandex 0.111 0.026 0.022 0.131
Return in-sample 0.055 0.227 0.172 0.055
Risk in-sample 0.245 0.218 0.245 0.218
Sharpe ratio in-sample 0.223 1.042 0.701 0.250
Return in 2017.07–2018.06 0.202 0.164 0.223 0.164
Risk in 2017.07–2018.06 0.179 0.170 0.215 0.153
Sharpe ratio in 2017.07–2018.06 1.129 0.965 1.039 1.067

This is how these four Russian portfolios look out of sample.

Returns for different optimal Russian portfolios out-of-sample

If one had invested in these four portfolios a year before, this is what one would have harvested by now.

Cumulative growth of one's investments into Russian stocks

So actually if one used this investing strategy based on the information available in mid-2017, by now, depending on the criterion (Sharpe ratio or maximum profitability) they would have earned 14–22% on Russian stocks or 15–17% on American stocks!

TODO: constrained optimisation using the best stochastic solution as a starting point.

<...>

Why are you still here? Why are you not trading yet?

Почему так важна нормальность ошибок в регрессии

В эконометрике очень полезно, если ошибки имеют нормальное распределение. В таком случае эконометристы сразу могут получить надёжные оценки для ошибок коэффициентов и интервалы для коэффициентов. При этом студентов обычно коварно успокаивают, мол, «если у вас очень много наблюдений, то тогда нормальность появится сама». Хо-хо, чёрта с два. Если у меня величина имеет \(t\)-распределение или, того хуже, равномерное, то нормальности у меня автоматически никакой не появится. Тогда гегемоны-семинаристы говорят: «Ой, ну ладно вам придираться, был доверительный интервал шириной 2,8 — стал шириной 2,9».

Сегодня я покажу, что отсутствие нормальности приводит к куда более страшным последствиям.

Рассмотрим функцию Кобба—Дугласа, и пусть, как всегда, отдача от труда равна \(\frac13\), от капитала — \(\frac23\), производственная функция описывается формулой
\[Y = A \cdot L^\alpha \cdot K^\beta \cdot f(\varepsilon). \]

Очень важно определить, что такое \(f(\varepsilon)\). Если мы хотим получить линейную модель вида
\[\ln Y =\ln A + \alpha \ln L + \beta \ln K +\varepsilon, \]
то при потенцировании должно быть
\[Y = A \cdot L^\alpha \cdot K^\beta \cdot e^{\varepsilon}.\]

Если \(\varepsilon\sim \mathcal{N}(0; 1)\), то как выглядит распределение \(e^{\varepsilon}\)? Давайте посмотрим.

В gretl сгенерируем Add — Random Variable — Normal со средним 0 и стандартным отклонением 1 (команда randgen(N,0,1)). Пускай для надёжности это будет десять тысяч, сто тысяч, миллион случайных величин! Нам же всего лишь график посмотреть, правда?

Summary statistics, using the observations 1 - 1000000
for the variable 'epsilon' (1000000 valid observations)

  Mean                     -0.0013720
  Median                  -0.00026809
  Minimum                     -4.9740
  Maximum                      4.7117
  Standard deviation           1.0007
  C.V.                         729.32
  Skewness                 -0.0015528
  Ex. kurtosis              0.0075161
  5% percentile               -1.6493
  95% percentile               1.6438
  Interquartile range          1.3480
  Missing obs.                      0

Test for normality of epsilon:
 Doornik-Hansen test = 2.76293, with p-value 0.251211
 Shapiro-Wilk W = 0.999998, with p-value 0.999982
 Lilliefors test = 0.000685939, with p-value ~= 0.3
 Jarque-Bera test = 2.75571, with p-value 0.252118

Normal distribution

Normal distribution 2

А вот и экспонента: Add — Define new variable: expepsilon=exp(epsilon).

Summary statistics, using the observations 1 - 1000000
for the variable 'expepsilon' (1000000 valid observations)

  Mean                         1.6477
  Median                      0.99973
  Minimum                   0.0069156
  Maximum                      111.24
  Standard deviation           2.1661
  C.V.                         1.3146
  Skewness                     6.1000
  Ex. kurtosis                 91.825
  5% percentile               0.19218
  95% percentile               5.1750
  Interquartile range          1.4510
  Missing obs.                      0

Exponent of normal distribution

Exponent of normal distribution 2

Так, это хорошо. Теперь представим, что после долгих лет упорного дата-майнинга (data mining) мы достали данные о миллионе рабочих и миллионе рабочих мест, оборудованных капиталом. Теперь сгенерируем рабочую силу и капитал для функции Кобба—Дугласа. Пусть они все будут равномерны, труд — от 3 до 15 (randgen(u,3,15)), капитал — от 5 до 20 (randgen(u,5,20)). Ничего, регрессия только выиграет от такой независимости. Add — Random variable — Uniform. Почему такие значения? Во-первых, все логарифмы больше единицы. Во-вторых, чем шире разброс регрессоров, тем лучше. Корреляцию необходимо тоже проверить статистически.

Да, и не забываем менять seed! Tools — Seed for random numbers. А ещё лучше записывать, какие числа используем.

Summary statistics, using the observations 1 - 1000000
for the variable 'L' (1000000 valid observations)
  Mean                         8.9972
  Median                       8.9914
  Minimum                      3.0000
  Maximum                      15.000
  Standard deviation           3.4639
  C.V.                        0.38499
  Skewness                  0.0028195
  Ex. kurtosis                -1.2002
  5% percentile                3.6035
  95% percentile               14.402
  Interquartile range          5.9986
  Missing obs.                      0

Summary statistics, using the observations 1 - 1000000
for the variable 'K' (1000000 valid observations)
  Mean                         12.500
  Median                       12.501
  Minimum                      5.0000
  Maximum                      20.000
  Standard deviation           4.3319
  C.V.                        0.34655
  Skewness                4.5311e-006
  Ex. kurtosis                -1.2022
  5% percentile                5.7482
  95% percentile               19.247
  Interquartile range          7.5089
  Missing obs.                      0

corr(L, K) = -0.00050070
Under the null hypothesis of no correlation:
 t(999998) = -0.500703, with two-tailed p-value 0.6166

Uniform factor distribution 1

Uniform factor distribution 2

На диаграмме разброса покажем только первые десять тысяч точек, так как оно всё рябит.

Factor distribution

Выполнены почти все условия теоремы Гаусса—Маркова ((1) отсутствие корреляции ошибок благодаря качественному генератору случайных величин (статистика Дарбина—Уотсона), (2) детерминированность и независимость \(X_i\), (3) отсутствие систематической ошибки \((\mathbb{E}(\varepsilon)=0)\), (4) гомоскедастичность ошибок (в силу, опять же, генератора)) — осталось только заняться спецификацией.

Всё отлично. Теперь придумываем модель. Пусть
\[Y = 2{,}71828 \cdot L^{0{,}333} \cdot K^{0{,}667} \cdot e^{\varepsilon},\]
\[Z = 2{,}71828 \cdot L^{0{,}333} \cdot K^{0{,}667} + 5\varepsilon.\]

В данном случае мы немного чувствуем себя богами, и это оттого, что нам известны истинные \(\varepsilon\). Что же, генерируем желаемые переменные в Add — Define new variable:

Y=2.71828*L^0.333*K^0.667*expepsilon
Z=abs(2.71828*L^0.333*K^0.667+5*epsilon)

Почему я применил abs() в последней формуле? Потому что кое-где эти дурацкие нормальные ошибки уводят наш чудесный выпуск в минус. На миллион наблюдений и при моих параметрах модели таких случаев едва ли наберётся сто, поэтому можно успешно их увести в положительную область. Всё равно они по модулю маленькие — действительность отражают неплохо.

Смотреть там особо не на что, пока мы не получим результаты регрессионного анализа. Выделяем все четыре. Нажимаем Variables — Add — Logs of selected variables. Теперь идём в Model — Ordinary Least Squares и задаём параметры. Зависим l_Y, независимы l_K и l_L.

Model 1: OLS, using observations 1-1000000
Dependent variable: l_Y
             coefficient   std. error   t-ratio   p-value
  -------------------------------------------------------
  const       1.00204      0.00811912    123.4    0.0000  ***
  l_L         0.333701     0.00229338    145.5    0.0000  ***
  l_K         0.665010     0.00261911    253.9    0.0000  ***

Mean dependent var   3.340627   S.D. dependent var   1.042599
Sum squared resid     1001331   S.E. of regression   1.000667
R-squared            0.078822   Adjusted R-squared   0.078820
F(2, 999997)         42783.15   P-value(F)           0.000000
Log-likelihood       -1419603   Akaike criterion      2839213
Schwarz criterion     2839248   Hannan-Quinn          2839223
Log-likelihood for Y = -4.76023e+006

Всё отлично. Правда, R² низкий. Но это ещё ничего не значит. Как видите, этот тест отлично показал, что R² — это вообще не индикатор. Это регрессия начинает извиняться, что в каждом наблюдении у неё есть ошибка, а ошибок много (очень много), а они ещё все скачут экспоненциально (вернее, дают большую ошибку в оригинальный \(Y\) и оставляют этот гадкий след даже после логарифмирования). Коэффициенты, конечно же, точные. Дополнительно исчисленный Дарбин—Уотсон (в таблице не показан) равен 1,998 792. Теперь вторая модель.

Model 2: OLS, using observations 1-1000000
Dependent variable: l_Z
             coefficient   std. error    t-ratio   p-value
  --------------------------------------------------------
  const       0.860155     0.00174908     491.8    0.0000  ***
  l_L         0.350244     0.000494058    708.9    0.0000  ***
  l_K         0.700475     0.000564230   1241      0.0000  ***

Mean dependent var   3.320815   S.D. dependent var   0.375988
Sum squared resid    46470.85   S.E. of regression   0.215571
R-squared            0.671274   Adjusted R-squared   0.671273
F(2, 999997)          1021021   P-value(F)           0.000000
Log-likelihood       115526.5   Akaike criterion    -231047.0
Schwarz criterion   -231011.5   Hannan-Quinn        -231037.2
Log-likelihood for Z = -3.20529e+006

Что?! Что я вижу?! На миллионе наблюдений в допущении экспоненты нормальных ошибок мы получили асимптотически смещённые оценки! Позор! Позор! Позор этому миру. Можно уже начинать сжигать книги.

Проведём тест на сумму коэффициентов. Я понимаю, что на миллионе в первой регрессии и так уже коэффициенты близки к истинным, поэтому не так хорошо видно, что́ я на самом деле пытаюсь понять, ибо \(0{,}333\,7+0{,}665\,0\approx\ldots1\) — а насколько именно приближённо? Вот на эту единицу мы и проверим. Tests — Linear restrictions. Вводим мантру:

b[l_L]+b[l_K]=1

Ответ для регрессии \(Y\):

Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 999997) = 0.137098, with p-value = 0.711183

Restricted estimates:
             coefficient   std. error   t-ratio   p-value
  -------------------------------------------------------
  const       0.999065     0.00116651    856.5    0.0000  ***
  l_L         0.334261     0.00172455    193.8    0.0000  ***
  l_K         0.665739     0.00172455    386.0    0.0000  ***
  Standard error of the regression = 1.00067

Смотрим p-value: 0,711 — конечно, при миллионе-то наблюдений! Отлично, тест пройден.

А если ввести линейные ограничения в модели с кривыми ошибками (регрессия \(Z\))?

Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 999997) = 4569, with p-value = 0

Restricted estimates:
             coefficient   std. error    t-ratio   p-value
  --------------------------------------------------------
  const       0.977157     0.000251871   3880      0.0000  ***
  l_L         0.328230     0.000372363    881.5    0.0000  ***
  l_K         0.671770     0.000372363   1804      0.0000  ***
  Standard error of the regression = 0.216063

Удивительно, но наши смещённые коэффициенты чудесным образом «пофиксились» самостоятельно! Другое дело, что p-value равно нулю, что мы улетели в «хвост» распределения, гипотеза о равенстве суммы единице отвергается (ошибка первого рода), поэтому данный тест показывает нам жирную дулю, что ненормальные ошибки — бич регрессий.

Давайте сохраним остатки (Save — Residuals) и посмотрим на их распределение. Заодно и потестируем. Я назвал их eY и eZ соответственно.

Test for normality of eY:
 Doornik-Hansen test = 2.76177, with p-value 0.251356
 Shapiro-Wilk W = 0.999998, with p-value 0.999989
 Lilliefors test = 0.000695252, with p-value ~= 0.28
 Jarque-Bera test = 2.75455, with p-value 0.252266

Good residuals 1

Good residuals 2

Test for normality of eZ:
 Doornik-Hansen test = 994212, with p-value 0
 Shapiro-Wilk W = 0.913314, with p-value 1.50967e-149
 Lilliefors test = 0.0623541, with p-value ~= 0
 Jarque-Bera test = 2.59783e+007, with p-value 0

Bad residuals 1

Bad residuals 2

Очень похоже на распределение Стьюдента, но я не могу знать наверняка, так как и пик острее и выше, и «хвосты» тяжелее.

Для полного счастья запишем, как будут выглядеть истинные ошибки в модели с функцией Кобба—Дугласа и прибавленной ошибкой:
\[Y = AL^{\alpha}K^{\beta} + \varepsilon = AL^{\alpha}K^{\beta}\left(1 + \frac{\varepsilon}{AL^{\alpha}K^{\beta}} \right)\]

Так как последняя скобка довольно близка к единице (уж слишком на многое делят несчастный эпсилон), то действует выражение эквивалентности \(x \stackrel{x\rightarrow0}{\sim} \ln (1+x)\). Но я есмь человек, любящий точность, а оттого посмотрю на оба распределения. Как видно ниже, они почти сливаются.

error1=epsilon/(2.71828*L^0.333*K^0.667)
error2=ln(1+epsilon/(2.71828*L^0.333*K^0.667))

error1
error2
Bad error

Вообще говоря, ничто не может быть нагляднее, чем график. Особенно трёхмерный. Конечно, я не могу построить график с миллионом точек. Даже сохранить данные не могу (программа зависла, пришлось к данному моменту всё заново делать... чёрт). Поэтому построим график с уравнениями конечных регрессий по выборке размером в две тысячи наблюдений.

Для «правильной» регрессии l_Y график выглядит так:
Volumetric normal distribution 1

Volumetric normal distribution 2

Volumetric normal distribution 3

Видно, что «пух» есть, что сильный, что на миллионе наблюдений он вообще превращается в персидский ковёр, но... Жизнь есть жизнь.

Для регрессии с неправильно специфицированными ошибками график выглядит так:

Volumetric not normal distribution

Volumetric not normal distribution 2

Volumetric not normal distribution 3

«Пуха» меньше, но и цена этого выше.

Вывод. Допущение о «нормальности» этого мира — абсурд. Несмотря на «нормальность» распределения ростов, на стьюдентовость распределения масс и на то, что уже на двадцати пяти монетках можно получить нормальное распределение, глупо, смешно и опасно предполагать, что парочка безобидных преобразований просто увеличит стандартные ошибки. Даже если это утверждал сам Андрей Александрович Мамонтов, посмеиваясь над скептически настроенными студентами. Правильная спецификация регрессоров даже асимптотически не может исцелить от ненормальности ошибок. Memento mori.

Постскриптум. «Иногда малая выборка лучше большой» — это верно или неверно? Неверно, но доказательство очень трудное. Тогда я скажу по-другому. Можно ли подпортить картину торжества модели с правильной спецификацией? Можно. Построим при помощи МНК те же две регрессии, но только на выборке из 800 наблюдений. Так как выборка — случайная величина, то и бояться её нечего.

Model 1: OLS, using observations 1-800
Dependent variable: l_Y
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       0.755572     0.294980      2.561    0.0106    **
  l_L         0.309048     0.0845139     3.657    0.0003    ***
  l_K         0.796715     0.0943801     8.442    1.47e-016 ***

Mean dependent var   3.375052   S.D. dependent var   1.057215
Sum squared resid    806.8346   S.E. of regression   1.006151
R-squared            0.096536   Adjusted R-squared   0.094269
F(2, 797)            42.58025   P-value(F)           2.69e-18
Log-likelihood      -1138.554   Akaike criterion     2283.107
Schwarz criterion    2297.161   Hannan-Quinn         2288.506
Log-likelihood for Y = -3838.6

Model 2: OLS, using observations 1-800
Dependent variable: l_Z
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       0.786741     0.0624576     12.60    2.61e-033 ***
  l_L         0.348992     0.0178945     19.50    1.48e-069 ***
  l_K         0.732322     0.0199836     36.65    4.14e-173 ***

Mean dependent var   3.333156   S.D. dependent var   0.379279
Sum squared resid    36.17173   S.E. of regression   0.213037
R-squared            0.685294   Adjusted R-squared   0.684505
F(2, 797)            867.7622   P-value(F)           8.2e-201
Log-likelihood       103.3827   Akaike criterion    -200.7654
Schwarz criterion   -186.7115   Hannan-Quinn        -195.3665
Log-likelihood for Z = -2563.14

Какая регрессия кажется лучше? Коэффициенты одинаково кривоваты и почти равны (и это при 800 наблюдениях! Вы там думали, что склеить пять наборов по 24 точки — это асимптотика, но нет, взгляните горькой правде в глаза), но у второй куда больше дьявольский R², да и ошибки коэффициентов поменьше... Однако отсутствие асимптотики не мешает второй модели по-прежнему запарывать тест линейного ограничения:

(Регрессия Y)
Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 797) = 0.702842, with p-value = 0.402082

Restricted estimates:
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       1.00047      0.0410358    24.38     1.49e-098 ***
  l_L         0.262031     0.0632130     4.145    3.76e-05  ***
  l_K         0.737969     0.0632130    11.67     3.49e-029 ***
  Standard error of the regression = 1.00596

(Регрессия Z)
Restriction:
 b[l_L] + b[l_K] = 1
Test statistic: F(1, 797) = 9.2669, with p-value = 0.00241004

Restricted estimates:
             coefficient   std. error   t-ratio    p-value 
  ---------------------------------------------------------
  const       0.975023     0.00873522   111.6     0.0000    ***
  l_L         0.312843     0.0134561     23.25    1.05e-091 ***
  l_K         0.687157     0.0134561     51.07    1.13e-253 ***
  Standard error of the regression = 0.214138

Постсубскриптум. Вы только что прослушали «Фантазию» для двух фортепиано А. Н. Скрябина не сравнение двух моделей, а сравнение двух спецификаций, так как в обоих случаях регрессоры были одними и теми же.

Забавный фрагмент лекции Арефьева

Экономический рост. Читает Николай Геннадьевич Арефьев, записывает и дополняет Ваш покорный слуга.

Загадка премии за образование. В США и Европе были изменения в 1970-е, сейчас они происходят в России. Повысилась доля людей, получающих высшее образование. В 1950-х годах в США высшее образование получали 30 %, в 90-х — 75 %. В России сейчас в когорте 20-летних доля высокообразованных увеличилась. Значит, надбавка за образование должна уменьшиться. Но она выросла. И этому есть много объяснений.

Эффект размера рынка. Пусть есть 2 блага. \(x_s\) — skilled, \(x_u\) — unskilled. Два человека. Вышкинец и птушник.

\[ \displaystyle X_s=\int\limits_0^1 A_{is} x_{is}^\alpha\,di \cdot L_s^{1-\alpha},\quad X_u=\int\limits_0^1 A_{iu} x_{iu}^\alpha \, di \cdot L_u^{1-\alpha} \]

\(x_{iu}^\alpha\) с коэффициентом 0,1 — это лопата, с коэффициентом 0,2 — рукавицы, 0,3 — водка, 0,4 — сапоги для гуляния по грязи и плаванья в навозе. Интеграл — это знак суммирования. Монополист, производящий \(x\), может произвести только одну единицу, но без издержек. \(x_{is}=1\), \(x_{iu}=1\). С такими упрощениями можно переписать производственную функцию.

\[\displaystyle X_s=A_s L_s^{1-\alpha},\quad X_u=A_u L_u^{1-\alpha},\quad A_s=\int\limits_0^1 A_{iu} \, di\]

\[\displaystyle w_s=p_s (1-\alpha) A_s L_s^{-\alpha},\quad w_u=p_u (1-\alpha) A_u L_u^{-\alpha}\]

Что есть премия за эдьюкейшн? Сие есть во сколько раз зарплата умного вышкинца выше зарплаты тупого птушника.

\[\displaystyle \frac{w_s}{w_u} =\frac{p_s}{p_u},\quad \frac{A_s}{A_u} = \frac{L_s^{-\alpha}}{L_u^{-\alpha}}\]

Знакомы с функцией постоянной эластичности замещения? Нет?! Facepalm. Это CES: \(u=(x_1^\beta+x_2^\beta )^{\frac{1}{\beta}}\), \(\beta<1\). Причём \(\beta\) отрицательным быть может. При \(\beta\rightarrow0\) \(u\rightarrow\sqrt{x_1 x_2}\). При \(\beta\rightarrow 1\) получаем совершенные субституты. При \(\beta\rightarrow-\infty\) получаем совершенные комплементы.

Толстовки с логотипом ВШЭ стали известны благодаря письму

Для меня стало приятной новостью известие о популярности моего электронного сообщения производителям толстовок с логотипами вузов, в том числе и ВШЭ. Поэтому привожу далее оригинальный текст без изменений.

//////////////////////////////
Привет Эконому!

Мы решили, что нашему факультету не хватает одежды, которая говорила о нашей принадлежности к главному факультету НИУ-ВШЭ. Поэтому мы решили организовать закупку толстовок и джемперов с надписью "ECONOM" и логотипом ВШЭ. Мы сделали их качественными и недорогими. (http://vk.com/tolstovki_econom)Доступные цвета : розовый, красный, серый и белый.Надпись и логотип сделаны вышивкой высочайшего качестваТкань толстовки/джемпера очень качественная, при стирке не образуется катышков. Зимой не будет холодно, а летом жарко.

Девчонки, толстовки не будут висеть на вас=) имеются как мужские, так иь женские размеры=)

Цена толстовки 1380р

Цена джемпера 1220р

!!! Толстовки можно посмотреть и примерить в Вышке. 

При желании можно заказать толстовку/джемпер без вышивки логотипа и надписи. (860р за толстовку и 720р за джемпер) Значительно дешевле, чем в магазине=)

По всем вопросам +79165494140 Артур

Также адрес странички с более полной информацией http://vk.com/tolstovki_econom

Уважаемый Артур,

Да, факультету действительно не хватает одежды. Особенно летом. Только наступает июнь — почему-то внезапно все перестают носить толстовки, и потные небритые ляжки начинают ужасать взор даже тех, кто экономистом не является (что уж говорить о тщедушных апологетах ойкоса!). Рахитичные и далеко не модельной внешности студенты своею наготой дискредитируют факультет, а им, тёмным, никто даже ничего не подскажет из-за заискивающей толерастии.

Спасибо, что уделили внимание текстильному вопросу и охарактеризовали товар. Но всё-таки моё мнение диаметрально противоположно Вашему. Позвольте оное аргументировать.

Я уже однажды отзывался в письменной форме о данных толстовках, но разрешите ещё раз изложить эти тезисы сжато. Качество шрифтовой работы хромает в первую очередь. Об оптическом кернинге знаков речи даже не идёт. Студенты ВШЭ выгодно отличаются от прочей массы студентоты именно тем, что способны примечать малейшие изъяны и их конструктивно критиковать. Настоящий экономист усмотрит в применении простого дугообразного (arc) преобразования траектории текста баскетбольной гарнитурой из Adobe Illustrator халтуру и безыскусность. Он не побоится об этом сказать.

Подобная порочная практика «маркировки» студента — это выявление его предпочтений. Пусть гарвардские белые воротнички и будущие нервные клерки (минус очко тем, кто прочитал это как «нервные клетки») напяливают на себя клеймо — настоящий экономист никогда не предоставит информации о своей системе предпочтений производителю во избежание дискриминации; если студенты Гарварда избирают путь самонасилия и рабства, то это разруха в их головах. Культ факультета или вуза — это то же самое, что и культ бренда. Демонстративное потребление, отжимающее у консьюмтариата то, что делает его заработную плату эффективной, обезличивает и смешивает массы. Кто доволен беззащитностью жертв компании Apple перед любой нестандартной задачей, не предусмотренной разработчиками? А кто в оной виноват? Сами служители культа.
Те, кто будет проходить по улице мимо носителя данной толстовки из очень качественной ткани, будут произносить одну из трёх фраз:
— Снова размундирились, всё какие-то однодневные словечки зарубежные да символы загогульные носят на фуфайках! (Не знающие о ВШЭ.)
— Опять эта раздутая ВШЭ, экономистов как собак нерезаных, рекламирует себя — аж трещит, уже и молодёжь превратила в доски объявлений! (Просто знающие о ВШЭ.)
— Ему во втором модуле было мало, так он ещё и в свободное время таскается с этой зубастой монограммой, сдобренной дилетантской типографикой! (Студенты ВШЭ.)
Три группы произносящих это составляют генеральную совокупность населения, поэтому ничего, кроме истого негодования, носители данных фуфаек не словят на недружелюбных стогнах града Москова.

Прочли ли дизайнеры данных толстовок хотя бы одну книгу о дизайне? О типографском наборе и шрифтовом искусстве? «Каллиграфия для всех»? «Краткие сведения по типографскому делу»? Жаль, ведь экономистам нравится чувствовать себя сэрами и аристократами, они бы пообщались с Вами на данную тему. Это лишний раз подчёркивает то, что экономист — это боевая машина, теург, спецназ, вихрь природы, звёздный ажиотаж и вся Вселенная. Почему на фуфайке не видно, что эконом — это центр Вселенной?

Логотипа ВШЭ не обругивал только ленивый. Почему вы идёте против инициативности? Почему бы не сотворить собственный логотип Вышки? В качестве бонуса присылаю скорый набросок. Если экономист — индивидуальность (а так оно зачастую и есть; правда, эти фрукты толстовок не носят, но об этом позволю себе скромно умолчать), то он не должен кряхтеть под гнётом ржавых стереотипов и готовых предрассудочных полуфабрикатов.

Higher School of Economics logo

И напоследок десерт из лексики, логики и грамматики.

  • В первом предложении пропущена частица «бы».
  • Кавычки в русском языке принято делать «такими (вложенные — „такими“)», а не "такими". Любая типографская раскладка для клавиатуры этот порок лечит (минус очко тем, кто прочитал это как «понос лечит») всего за два нажатия, а также позволяет употреблять в предложении корректный знак для тире — такой — а не мелкий-мелкий дефис (почувствовали разницу?).
  • Не вы сделали толстовки качественными и недорогими, а производственная технология, эффективная вертикальная структура и благосклонная к конечному потребителю схема ценообразования (надеюсь!), поэтому предложение надо переписать.
  • Перед двоеточием пробел поставлен некорректно, его быть не должно: таково правило.
  • Если в предложении употребляется смайлик, то его принято отделять от текста предшествующим пробелом. После него сразу ставится знак препинания, завершающий предложение. Смайлик не заменяет, а дополняет пунктуацию. Пример: ой ли :-)? Но лучше от использования смайликов отказаться, ведь мы не в зверинце живём, верно? Если автор хочет, чтобы адресат улыбнулся, то он должен не императивно его тыкать в командную комбинацию знаков («Смейся, адресатишка, здеся»), а искрить остроумием и литературной игрой. Смайлики происходят от неумения управлять родным языком, от нежелания научиться пробуждать в читателях определённые эмоции при помощи богатства словесных выражений.
  • Запись о цене должна выглядеть так:
    Цена толстовки — 1 000 руб.
    Если у кого-то система новее, чем Windows XP, то корректно отобразится более уместный знак тонкого пробела. (Цена толстовки — 1 000 руб.) Каждые три разряда разделаются пробелом или тонким пробелом, в качестве символа валюты употребляется установленное ЦБ символьное обозначение «руб.». Экономисты об этом знают.
  • Номер телефона для облегчения восприятия и разгрузки зрения принято записывать так: +7 916 549-41-40.

Уважаемые Артур и прочие предприниматели, если вы пришлёте в ответ полностью исправленное, написанное в соответствии с нормами типографского набора и русской письменной речи коммерческое предложение (основные допущенные вами ошибки описаны выше), то я куплю у вас толстовку.

С наилучшими пожеланиями,
Андрей Викторович Костырка — он не боится это сказать.

P.S. Эй, вы, моё предложение ещё в силе, слово твёрдо! Что же вы не отвечаете? Не хотите, чтобы я у вас фуфайку купил?

Что придёт в будущем на смену деньгам

Необходимо отказаться от всех мировых валют. Возврат к золотому стандарту лучше, но его историческая несостоятельность отпугивает нас от его введения.

Необходимо исчислить всё мировое богатство. В каких единицах? Ни в каких. Все богатства мира — это единица (1). Следующий шаг — поделить мировое богатство по материальным благам так, чтобы сумма их стоимостей давала единицу. Скажем, буханка хлеба — это 5,75 × 10−19 МБ (мирового богатства), мой настольный компьютер — 3,6 × 10−16 МБ. Не будет инфляции. Не будет валютных кризисов. Рост мирового богатства (скажем, прибавление стоимости технологических разработок) будет выражаться в относительном удешевлении того же хлеба. При этом соотношение цены хлеба и компьютера будет одинаковым что до, что после оформления дорогостоящей научной разработки.

Banks will go bankrupt! World government will go bankrupt! Problems, Friedman?

Исчезнут бумажки, правительства не смогут накачивать свою валюту с учётом рыночной конъюнктуры. Установится совершенная конкуренция в сфере обращения.

Запись будет отредактирована, расширена, превратится в глубокое научное изыскание.

Savoir-vivre от Адама Смита

Обожаю старинных писателей. Тогда то, что сейчас пишется для эпатажа, было зачастую обыденной вещью. Сейчас за такие слова писателя сочли бы фриком, извращенцем или модернистом, а тогда это писалось искренне, безо всякого предосуждения.

Бедность низших слоёв народа в Китае далеко превосходит бедность самых нищенских наций Европы. В окрестностях Кантона многие сотни, как обычно утверждают, даже тысячи семейств не имеют совсем никакого жилища на суше и живут постоянно в маленьких рыбачьих лодках по рекам и каналам. Пропитание, которое они добывают себе здесь, настолько скудно, что они жадно выуживают самые негодные отбросы, выкидываемые за борт европейских судов. Любая падаль, например дохлая собака или кошка, хотя бы совсем разложившаяся и испускающая зловоние, столь же лакомая пища для них, как самая здоровая пища для народа других стран. Браки поощряются в Китае не выгодой, получаемой от детей, а дозволением умерщвлять их. Во всех больших городах каждую ночь много детей оставляют на улице или топят, как щенят, в реке. Утверждают даже, что выполнение этого ужасного дела является признанной профессией, которая даёт пропитание многим людям.

Грека и европейский кризис

Ехал Грека через рынок,
Видит Грека — в рынке риск;
Сунул Грека риск свой в рынок
И разрушил рынок вдрызг!

Почему супруги должны консолидировать усилия на кухне

Последнее время слышится всё больше призывов к домохозяйкам: «Хватит быть рабыней!», «Перестаньте готовить мужу, пусть готовит сам!» и прочая. Мало ли на свете глупцов, несведущих в математической науке!

Польза приготовления человеком еды более чем на себя доказывается на конкретном примере. Производственная функция готовящего на кухне имеет свойство возрастающей отдачи от масштаба, т. е. \(f(kx,ky)>kf(x,y)\). Это означает, что, для того чтобы приготовить в два раза больше пищи, не нужно затрачивать вдвое больше времени (или иного ресурса). Поэтому с целью ублажения сумасшедших и безголовых феминисток я предложу вариант, который устроит всех без исключения: пусть жена с мужем готовят еду на обоих по очереди, так как, скажем, если на приготовление одной порции салата или супа уходит \(t\) единиц усилий и \(r\) единиц времени, то на приготовление двух порций придётся потратить лишь около \(1{,}8t\) усилий и \(1{,}5r\) времени. Это связано с тем, что утварь, кухонные приборы и прочие средства производства будут обслуживаться лишь один раз (достать, помыть), и поэтому для обработки двойного объёма ресурсов требуется ещё меньше времени, ведь суп может вариться в ёмкости большего объёма.

Имеем: при полном разделении труда и индивидуальности членов семьи и муж, и жена в сумме тратят \(2t+2r\) усилий и времени, в то время как вступившие во взаимодействие благоразумные их соседи — лишь \(1{,}8t+1{,}5r\), что позволяет им сэкономленные силы и время перераспределить с учётом переговорной силы сторон.

Где ваш рассудок, борчихи за права женщин?