text:
# What was Burroni's sketch for topological spaces?

In a 1981 talk, René Guitart cites Albert Burroni as having given "A first interesting example of a mixed sketch...for the category of topological spaces" in 1970. This was apparently done in Burroni's Université de Paris thesis, and I don't see any evidence of it having ever been published more, well, publically.

It seems that this must be a large sketch, since categories sketchable by a small sketch are accessible, and in any case topological spaces would be locally presentable if they were accessible, and thus sketchable by a limit sketch.

Regarding large sketches, one has the large limit sketch of Spanier giving quasi-topological spaces. Burroni's thesis title includes the word "quasi-topologies," although it's not clear to me whether that is meant in Spanier's sense. Could Burroni have done something similar, adding some colimit cones to Spanier's topology on compact Hausdorff spaces to restrict to precisely ordinary spaces?

In any case, does anyone know his construction, or another sketch for spaces, or a reference in which one is written down? The only comparable thing I can think of is (EDIT:) Barr's theorem on spaces as relational $\beta$-algebras, ($\beta$ being the ultrafilter monad) and I doubt there's going to be some theorem relating relational algebras and large-sketchable categories.

• Actually, the relational $\beta$-module characterization is due to Barr. It is of course related to the Manes characterization of compact Hausdorff spaces as (functional) $\beta$-algebras. – Todd Trimble Jul 26 '17 at 21:28

If you read french, you can look at page 6 here. It is far from formal but it gives a good idea of the mixed sketch. In my understanding, there is a unique inductive cone in that mixed sketch, which comes from the canonical representation of the filter endofunctor as a colimit of representables. If you read on the rest of the paper, you'll see that this is a quite powerful method as Robert Paré showed that it could serve as a guide to sketch any category with split idempotent.

Here is a rough translation:

When trying to sketch topologies, the difficulties boil down to the problem of getting continuous applications as the natural transformations between models of the sketch. This prevents us to define a topology from its set of open sets, as it would eventually yields open maps instead of continuous one. Inspired by the idea behind Kowalsky's quasi-topologies, for which Ehresmann showed some interest, Albert [Burroni] proves that a topology on a set $E$ is no more than a "convergence" relation [denoted "$\to$"] between the set $\mathbf F(E)$ of filters on $E$ and the set $E$ itself, required to satisfy some axioms. A continuous map is then a map that preserves this notion of "convergence". Among many equivalent others, here is a system of axioms Albert proposed. The identifiers match the ones he used:

• $\mathbf T_2$ if $F \to x$ and $F\subset F'$, then $F' \to x$
• $\mathbf T_1$ $I_E(x) \to x$ ($I_E(x)$ stands for the principal filter on $x$)
• $\mathbf T'_4$ if $\phi\to F$ and $F\to x$, then $K_E(\phi) \to x$ ($K_E(\phi)$ stands for the filter on $E$ obtained by "summing" the elements of the filter $\phi$ on $\mathbf F (E)$)

By describing abstractly this process, he obtains the wanted sketch. It is formed by an "object of points" $e$, an "object of filters" $f(e)$, and an "object of filters on filters" $f(f(e))$, arrows from $e$ to $f(e)$ and from $f(f(e))$ to $f(e)$ to represent the maps $I_E$ and $K_E$, and a monomorphism with target $f(e)\times e$ to represent the convergence relation between filters and points. It is straightforward to translate the axioms above, after which it is just a matter of "typifying" the object $f(e)$ as the set of filters on $e$. To do so, Albert has the idea to decompose the foncteur $\mathbf F : \mathsf{Set} \to \mathsf{Set}$ as a canonical inductive limit of representable functors and then adds the necessary projective (discrete) cones and inductive cone to the skecth. The set-theoritical model of the sketch obtained this way are precisely the topological spaces. Notice that Albert also studies the model of this sketch in $\mathsf{Top}$ to conclude, apparently with a bit of disappointment:

"A topological topology is nothing else than a couple $(\pi,\pi')$ of topologies on a shared set $\theta(\pi)=\theta(\pi')$." [Burroni's thesis, p.65]

• Ok, great, it's essentially given as a relational algebra, for the functor of filters instead of ultrafilters. Thanks very much for the reference! – Kevin Carlson Aug 4 '17 at 16:09
• Possibly a dumb question, but what's the domain and codomain of the filter endofunctor and how does it act on morphisms? – goblin Feb 24 '19 at 23:59

The category of topological spaces is the category of models of a relational universal strict Horn theory $T$ without equality, i.e. the axioms are of the form $(\forall x)(\phi(x) \rightarrow \psi(x))$ where $x$ is a set of variables and $\phi$ and $\psi$ are conjunctions of atomic formulas without equality. The result was proved for topological spaces first by E. Manes in Equational theory of ultrafilter convergence, Alg. Univ. 11 (1980), 163-172 and generalized by J. Rosický to any fibre-small topological category in Concrete categories and infinitary languages, J. Pure Appl. Alg 22 (1981), 309-339. And indeed, the class of axioms is a proper class, otherwise the category of topological spaces would be locally presentable, which is known to be false. Indeed, no non-discrete space is presentable.

• Thanks very much, Philippe! I suppose that this constructs spaces as the models of a non-small limit sketch. So I still wonder why Burroni's sketch was cited as an interesting example of a mixed sketch. – Kevin Carlson Jul 27 '17 at 4:42
• @KevinCarlson Yes this proves that the topological spaces are the models of a large (i.e. non-small) limit sketch. My answer was about your very last question. You can find an explicit description of the axioms in Example 5.4 of J. Rosický's paper (which is in open access, you just have to google the title). – Philippe Gaucher Jul 27 '17 at 7:47
• Ah, great, I understand. – Kevin Carlson Jul 27 '17 at 16:59

The reference for Burroni's sketch of topological spaces is Comptes rendus de l'académie des sciences, 27 juillet 1970 (tome 271), which amazingly can be found online.

Perhaps Albert's thesis, Esquisses des catégories à limites et des quasi-topologies is still available. If you go to this Cahiers page and follow the link on the left hand page to ESQUISSES MATHEMATIQUES, it says

Most of them are still available and can be obtained by sending an e-mail

to the email address given there. Albert's thesis is no. 5 on the list.


text:
## Calculus: Early Transcendentals 8th Edition

a) 69.6 $\frac{heartbeats}{minute}$ b) 71.75 $\frac{heartbeats}{minute}$ c) 71 $\frac{heartbeats}{minute}$ d) 66 $\frac{heartbeats}{minute}$ Conclusion: We can find the slope at t = 42 by the average of the closest slopes which are t = 40 and t = 44, so m = $\frac{71+66}{2}$ = 68.5 $\frac{heartbeats}{minute}$
To find the slopes calculate by: m = $\frac{\frac{heartbeats}{minute}_{tf} -\frac{heartbeats}{minute}_{ti} }{tf-ti}$ a) m = $\frac{2948-2530}{42-36}$ b) m = $\frac{2948-2661}{42-38}$ c) m = $\frac{2948-2806}{42-40}$ d) m = $\frac{3080 - 2948}{44 - 42}$

format: plain


text:
cardinality constrained mean reversion - proof of concept

Here is an attempt to implement the ideas in this research paper titled "Identifying small mean reverting portfolios": http://www.di.ens.fr/~aspremon/PDF/MeanRevVec.pdf

The gist of the paper is that given N stocks in a universe, the algorithm identifies k stocks that are mean reverting (stationary). As a proof of concept, the attached code does not implement a trading strategy but only logs if the smaller portfolio is stationary or not.

Unfortunately, the optimization using cvxopt for this algo is quite intensive and quantopian blows up with out of memory or timeouts if universe contains more than 15 or so stocks. I think this algo will work well if universe is large (>50 stocks) and we select a basket of 8-10 stocks that are mean reverting.

Please play around with it and let me know if you identify any defects.

168
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
# Backtest ID: 553a868cd40f100d9417a1e3
There was a runtime error.
49 responses

A simple crude trading strategy around the stationary series.

168
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
# Backtest ID: 553a9323ef02fe0d1a37a387
There was a runtime error.

Satyapravin, This is a really interesting result. Looks like all the 'jumps' in the portfolio are to the upside. Have you tried backtesting it over a longer time frame or over other stocks or ETF's?

Disclaimer

The material on this website is provided for informational purposes only and does not constitute an offer to sell, a solicitation to buy, or a recommendation or endorsement for any security or strategy, nor does it constitute an offer to provide investment advisory services by Quantopian. In addition, the material offers no opinion with respect to the suitability of any security or specific investment. No information contained herein should be regarded as a suggestion to engage in or refrain from any investment-related course of action as none of Quantopian nor any of its affiliates is undertaking to provide investment advice, act as an adviser to any plan or entity subject to the Employee Retirement Income Security Act of 1974, as amended, individual retirement account or individual retirement annuity, or give advice in a fiduciary capacity with respect to the materials presented herein. If you are an individual retirement or other investor, contact your financial advisor or other fiduciary unrelated to Quantopian about whether any given investment idea, strategy, product or service described herein may be appropriate for your circumstances. All investments involve risk, including loss of principal. Quantopian makes no guarantees as to the accuracy or completeness of the views expressed in the website. The views are subject to change, and may have become unreliable for various reasons, including changes in market conditions or economic circumstances.

The algo is too computationally intensive and demands lot of memory and CPU. I don't think it scales well in quantopian unfortunately. Have to run it outside quantopian with a larger universe of atleast 50-100 stocks.

There's an outfit (also in Boston) that offers high-performance computing for financial applications:

https://elsen.co/

Perhaps they could lend a hand, if the Q research platform can't handle it.

Thanks for the reference Grant. I was thinking of splitting the universe into buckets of 10-15 stocks and trying to identify a mean reverting subset in each basket on Q. I will take a look at elsen.

Finally managed to crack it. The lesson I learnt is that mean reversion is very noisy around mean and non gaussian. Standard deviation or zscores don't help. I need to find a better measure of dispersion. Anyway this backtest has a sharpe of 2 and can be improved because there is a long period of 6 months where algo does not trade at all. This can be changed by increasing the universe of stocks to scan for opportunities.

168
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
# Backtest ID: 553f0cdee0e8390d36933492
There was a runtime error.

Hi Satyapravin,

The paper http://www.di.ens.fr/~aspremon/PDF/MeanRevVec.pdf is quite involved; I'd really have to roll up my sleeves to understand it (although I have some of the basic linear alebra, statistics, and conceptual tools). Generally, my interest is to see if the approach could be applied in the Q research platform to identify baskets of securities that could then be used in an algorithm like OLMAR (http://arxiv.org/ftp/arxiv/papers/1206/1206.4626.pdf) or a variant.

I have the feeling that the Q research platform may "break down" at some point. For example, say I wanted to pluck out mean-reverting baskets of 50-200 stocks out of a universe of 2000. Would it work, if I'm willing to wait hours/days?

Would you be willing to provide a break-down of what portion of the paper you applied, what it does, and why it works?

Also, it seems important to be able to explore the time scale of mean reversion? Does this analysis provide a means to do that? For example, if the time scale is too short, then Quantopian, running on minute bars, won't work, and if it is too long, then the deviation will have to be very large (which seems unlikely).

Thanks,

Grant

Hi Grant,

The paper is complex but after several reads became much clearer. Here is the gist of the paper. First we run lasso regression on stock prices to obtain an estimator which is the coefficient in VAR(1) model like this:

for j in range(0, col):
Ahat[:,j] = Lasso(0.1).fit(pval[:-1,:], pval[1:,j]).coef_

Then we compute the covariance matrix using shrinkage model like this:

oas = OAS()
oas = oas.fit(pval)
B = oas.covariance_

The next step is to solve a semi definite optimization program described below. Most of the complexity is because Q does not support CVXPY and we have to make do with CVXOPT.

# minimize    Tr(MY)
# subject to
#           ||Y|| - kz <= 0
#           Tr(Y) - z   = 0
#           Tr(BY)      = 1
#           Y          >= 0 - positive semidefinite
#
# where
#           Y =    X             z =     1
#               --------             ----------
#                Tr(BX)                Tr(BX)
#
#           M = A'B A
#           A = lasso regression estimator
#           B = covariance matrix
#           X = xx'    - x are the weights
#           k = max no. of stocks in mean reverting portfolio
#

Important points:

1. Since X = xx' we require N square variables for N stocks in universe. This quickly explodes. For example if universe has 30 stocks then optimization requires 900 variables and Q dies.

2. Most times mean reversion quickly collapses out of sample.

3. I used adfuller to check for stationarity which implies mean reversion

4. Sometimes the weights are heavily skewed. What I mean is that if there are 10 stocks are we are investing 100K, then if a single stock has an exposure of 50K it is uneven. But that is easy to solve by introducing a new constraint on each stock's exposure.

5. You can use half life to measure the time it takes to mean revert. Please see fit_OU method.

6. There is also another paper which introduces variance in portfolio. That can be clubbed with this strategy by introducing a new constraint.

By the way, since mean reversion breaks out of sample, I am actually trying out a strategy that bets that mean reversion is going to break. That means if a stationary portfolio is 2 standard deviations about mean value then we assume that it will break out even more. This seems to work as well. I will post results shortly.

I think the algorithm in paper works well to identify stationary portfolios but I am not experienced to trade it. This is where I am failing.

Also it does not make sense to look for mean reversion in a portfolio of random stocks. A friend of mine who works in this industry suggested that I use clustering techniques to group similar stocks and then look for mean reversion amongst similar stocks. For example look for stocks that move in same direction with oil prices and build a mean reverting portfolio of these stocks.

Yes, brute force searching for mean reverting portfolios of random stocks is very time-consuming. On Quantopian, it's also really important to only use these techniques in really liquid assets with narrow spreads, otherwise you'll find phantom mean reversion due to bid-ask bounce in the trade data.

On that note, here is a clustering algorithm demo in python if it wasn't shared before. It uses affinity propogation (haven't understood it yet) on correlations to group stocks into clusters.

Thanks Satyapravin,

At least now when I refer back to the paper, maybe I can focus on the right things.

Regarding the clustering, I've attached an example that was posted on Quantopian some time ago. Perhaps you'll find it useful.

Grant

89
Backtest from to with initial capital
Total Returns
--
Alpha
--
Beta
--
Sharpe
--
Sortino
--
Max Drawdown
--
Benchmark Returns
--
Volatility
--
 Returns 1 Month 3 Month 6 Month 12 Month
 Alpha 1 Month 3 Month 6 Month 12 Month
 Beta 1 Month 3 Month 6 Month 12 Month
 Sharpe 1 Month 3 Month 6 Month 12 Month
 Sortino 1 Month 3 Month 6 Month 12 Month
 Volatility 1 Month 3 Month 6 Month 12 Month
 Max Drawdown 1 Month 3 Month 6 Month 12 Month
# Backtest ID: 55416fc8d8811d0d24656440
There was a runtime error.

@Satyapravin B., You appear to be a pretty smart guy. And I wish you all the best if you decide to participate in the Q's Open. But your code is so obscure as to be impenetrable to my tiny brain. I've reviewed a number of your submissions and frankly have to look away else my mind get sucked into some quantitative induced vortex from which I'll never escape. Sluuuuuurp -- there went my half a billion neurons. But I'm drawn to your investigations like a moth to a flame. Singed every time however, I wonder if there's not some salve to ease my pain. For example:

def getGs(n):
G = np.zeros((n**2,2*n**2+1))
row = 0
for i in range(0,n):
for j in range(0,n):
g = np.zeros((n,n))
if i == j:
g[i,j] = -1.0
else:
g[i,j] = -0.5
g[j,i] = -0.5
g = np.ravel(g).tolist()
g = g + ([0] *(n**2+1))
G[row,:] = g
row += 1
return [matrix(G)]

My head has started to throb just seeing this code from the periphery of my eyesight.

Is it at all possible to implore you to elaborate your code in such a way as to render it intelligible to the quantitatively challenged -- such as myself?

Hi Market Tech,

Apologies. I am trying hard to enter Q open but haven't managed to find a consistently profitable strategy. If things work in 2013-15 they fail miserably in 2006-2008 or something like that. I was in a rush to get something going and hence the lack of comments in my code.

Unfortunately CVXOPT is very complex. All of the complexity will reduce once Q supports CVXPY which I have requested. I will post a cleaner version as soon as its done.

Also I am not a python programmer which you would have known by now looking at my code. I program in C++ mostly and am learning python just for Q.

always great to be able to learn from someone like you Satyapravin, but I must say that it takes me a long time to read the lovely condensed code, lack of brain power like MarketTech I guess

@ Satyapravin

Regarding the contest, keep in mind that the only time periods that matter are the 2-year backtest, 1-month paper trading, and 6-month real-money trading. Personally, I didn't even bother backtesting beyond the 2-year period, and I got very lucky during the paper trading period, and so far I've been unlucky trying to profit off of Q's capital (which, by the way, I knew could happen--the contest is all about reward to the winners with all of the risk burdened by Q).

@ Market Tech,

As a professional programmer (right?), I think I hear you volunteering to assist Mr. Bezwada in writing elegant, computationally efficient code that could be read by a 5 year old. Regarding the example you provided, getGs(n), it appears to be a matrix of constants needed by the secret sauce optimizer, CVXOPT's sdp. If you really want to dispense with some neurons, have a look at http://cvxopt.org/userguide/coneprog.html?highlight=sdp#semidefinite-programming (or you could buy a nice bottle of Scotch instead). This is actually a great case for Q to consider in their quest for domination of the hedge fund space. It is my belief (which I've expressed to them) that they are missing a key element, which is that hard problems get solved by diverse multi-disciplinary teams, not by "geniuses" working in isolation. The person who can understand and explain CVXOPT's sdp probably won't be the same person who can implement it best in code. And you may need an experienced, creative, intuitive type to steer the whole effort in the right direction, but who would never be able to understand what's going on "under the hood." And someone who understands the nitty-gritty mechanics of trading, market microstructure, etc. And so on. I'm guessing that the best hedge funds have such teams working elbow-to-elbow.

Hmm, hard to choose, I do rather like scotch...

Collaboration is always an option that the Q has curiously provided through its web interface. I've only ever used it once.

I would mention, as I am wont to do, that the trading headspace is an odd duck. On the media extolled side we have "proprietary algos potentially worth millions of dollars" -- so why would you ever think of sharing them (collaborating)? On the other side, the reality side, we have, as you've explained Grant, single disciplined individuals who would benefit greatly by teaming up with others who could then, as a group, build much better systems and strategies. I've run into this dichotomy time and time again.

"I can't share my secret trading technique."
"Yeah, but you can't code it yourself so you'll have to share it to build it."
"But, but, someone will steal it!"
"Right, so keep it to yourself then so that it will never actually make anyone money."
"But, but..."

This philosophy may not be as prevalent here, it is after all a forum of algo writers. But I would point out how few of the leaderboard contestants actually participate here. So maybe it's more systemic than one might think. It's an intriguing situation.

But I would point out how few of the leaderboard contestants actually participate here. So maybe it's more systemic than one might think. It's an intriguing situation.

True. The contest did change things a bit.

But I would point out how few of the leaderboard contestants actually participate here. So maybe it's more systemic than one might think. It's an intriguing situation.

True. The contest did change things a bit.

Is there a more-or-less canned, tried-and-true brute force method for scoring random baskets of securities (from a limited universe) for their suitability to collectively mean revert? If so, it would be straightforward to write some code in Q research, just to see how it plays out. There is no parallel processing, and there are unknown memory limits, but it is my understanding that a notebook could run indefinitely just churning over data.

I suppose one could do the clustering first, followed by testing each cluster. Are there other divide-and-conquer techniques?

Any thoughts? Pointers?

This wanders some distance off Satya's path, but that's pretty much what the Directed Acyclical Graph DAG was intended to do, identify a group of securities that were highly likely to mean revert. Alternatively one could perform a simple "arrhg!, who's with me?" kind of selection by picking an anchor security and then locating every other security that was either correlated or cointegrated with it. That way you'd limit your search. Pick 10 or 20 consistent anchors and then run dickyfuller or leastsqrs on the universe against each one; picking the top N as your group. Could even be automated perhaps and included in a strat.

Last year, in MATLAB, I tested all baskets of 2, 3 and 4 ETFs out of a universe of ~200 ETFs I think on two years of daily trade data, using Johansen's method for cointegration, several different stationarity tests, and a quick backtest as a final filter. It ran on my Macbook Pro for about four months. So far, none of them have really held up to intraday backtesting with quote data, leading me to believe that most of the so-called mean-reversion was just bid-ask bounce noise.

Why does "bid-ask bounce noise" lead to false positives? And is there a way to filter it out?

I deal with it by back-testing these things using quote data, and making sure that my backtest algo buys the basket at the asks and sells at the bids. I am not sure how one could get rid of it when using trade data, you just have to restrict your universe to those stocks with very narrow spreads.

I drew a picture of how it is deceptive non-mean-reversion on this thread: https://www.quantopian.com/posts/inaccurate-max-drawdowns-on-no-trade-days-slash-periods

Well, I'll have to think about why the bid-ask spread would come into play. I guess you are saying that price swings would need to be outside the bid-ask range to be useful, but if you don't know the bid-ask, then how would you ever know signal from noise.

Yeah, well imagine you had a perfectly hedged basket of X and Y which was perfectly stationary when long +1 of one and short -1 of the other. Suppose the value of this spread is now, and forever, $5. Generally, it would awesome to buy the spread when it dipped below 5, and sell it when it poked above$5.

Now, suppose that X has a bid-ask spread of $1 and so does Y. That means that the spread of the basket is$2, and so the bid of the basket is $4 and the ask of the basket is$6. Now, through the random noise trading of X and Y, some trades will be at the ask and some trades will be at the bid. If you calculate your basket value from trade data alone, you'll see the value of your spread jumping from $4 to$6 and in between. A correctly optimized algorithm for trading this mean-reverting spread would be to buy it when it dips to $4, and sell it when it trades at$6.

However, this is totally not executable, unless you are market-making both legs, since it costs $6 to buy the spread and it costs$4 to sell the spread. All the motion is really just random noise.

Statistical security trading seems replete with exceptions to the rules of mean reversion trading. I've built and tested so many variations of pairs, triangular arbs and basket strategies that I'm not even sure why I bother anymore. For the record NONE of the clients we ever had traded securities in such a fashion. They either day-traded equities using highly proprietary strats, traded options using broader strats or, most often, didn't trade equities at all. 90% of our clients traded futures. And there they had much better success using mean reversion in pairs and small groups. Grain complex, metals, energies, all with multiple expirations have much more predictable behaviors.

if you don't know the bid-ask, then how would you ever know signal from noise.

Indeed!

I believe the issue lies in the fact that there is so little information in the price values securities. On the surface it seems tractable, analyzing historic OHLCV price to project the path of said price over the next few minutes or days. But the lack of true predictive information tells another tale. This along with the mechanics, as Simon points out, belie the truth of profitable executions.

How much information is really in a WIO (weighted index oscillator) that has come out of its tight standard deviation channel and jumped up a deviation or two? Oh, it's bound to revert, stock A will fall back down or stock B will lift to match its brother; so sell A and buy B. But then the relationship changes. WIO is now at a new level. The information we had was incomplete. Stock A signed a big contract or stock B got hit by an EPA lawsuit. Price just doesn't contain enough information to reliably, profitably trade.

That's where all these social data feeds and insider trading stats (the legal kind) and earnings and fundy and blah blah blah type information sources try to overlay price with something, anything that can give credence to the reason why stock A jumped while stock B fell. But all these extra information sources are so damn hard to leverage. I bet the Chicago Traders Convention will have a ton of clever characters talking about their strikingly new data sources.

@ Market Tech,

For a small number of securities, the risk will always be there. As an institution, it seems you are gonna want something like a portfolio of at least 100 securities, and a minimum of $100M in capital. Then, you can have$1M in each security (on average) to keep trading costs down, and all of the special events should average out. Of course, this assumes that one could find 100 securities for which there is enough inefficiency for there to be an overall arbitrage, right? Sounds kinda hard.

@ Grant,

Once I get access to research platform, I will try to run it locally on my PC to see if CVXOPT can work with 100 stocks. That should provide some insights.

Found this interesting paper on how to model mean reverting spreads. Since above algorithm already finds a stationary portfolio the next step would be to trade the spread. This paper describes an algorithm to model these spreads so that they can be traded. Still above me for the moment but hopefully in time I will be able to understand and verify it.

@ Simon,

Regarding this business of the trade prices bouncing randomly between the bid and ask, shouldn't this be evident in a histogram of prices? Wouldn't one see a bi-modal distribution, at least over short time scales? And for OHLCV minute bars, wouldn't the H & L values be be representative of the ask & bid prices, respectively? Or are you thinking more in the limit of thinly traded securities where perhaps bid-ask becomes kinda murky? I did some quick skimming, and I gather that if I, guy-off-the-street retail trader goes to buy, I'll pay the ask, and when I sell, I'll sell at the bid. So, I'm supposing that if I looked at something like SPY in the research platform, I'd see trades distributed around two values (e.g. something like two superimposed Gaussian distributions). Or no?

Not really, since the mid-point of the bid-ask spread moves throughout the day. It contributes to the appearance of mean-reversion on short timescales, and I suppose it affects the kurtosis of returns, but it won't be bi-modal afaik. The bid-ask is never really murky, it's always there, no matter how thinly traded the securities are.

Perhaps I'll make a research notebook that has a model...

Yeah, if the midpoint is noisy or drifts quickly, then the bid-ask effect in the reported trades will be all smeared out, I suppose. But if you factor out the drift then shouldn't you see a bi-modal trade price distribution for large enough spreads? Or maybe the spreads are always buried in the noise?

Plot a couple of WIOs and see if you think you can trade them. In the end what you end up trading is a single composite price-like-thing. You can fairly easily see these in TradingView... For instance, here's a few simple WIOs

You get the picture.

Now, those are daily prices of course. And most likely one would try to pairs/basket trade on a much tighter periodicity. But then you get into Simon's case of what's a spread and what's a bid/ask wiggle? So you'd have to back out your trading analysis time period to one where that wiggle gets swamped by the actual drift in the spread. And then what you've got to worry about is your WIO going directional. And ultimately that's what kills every pairs strat. The basket strats, as Grant points out, are where I've tended to concentrate my arb strats. But still, they have their directional tendencies too.

If you look at futures calendar spreads they seem to be more easily traded in this way:

(In a basket arb strategy one ends up realizing that the instruments you trade are effectively (n * n-1) / 2 WIOs (pairs strategies). 10 legs gives you 45 pairs that you wind up trading. Now, with that thought in mind you get the added information of the composite WIO, either a blend of 10 individual technical variations, or a grander more complex bundle of 45 WIOs. Trading 45 WIOs, cherry picking the widest aberrations among them, may be the best bet.)

Grant: I stand corrected, in my simulations, if the spread is much wider than the drift, there are detectable traces in the distribution of returns, but only if the spread size is constant. If the spread size is stochastic too, then it's completely smeared away...

@ Simon,

That's the picture in my mind. If the spread dominates, then it'll be detectable. But if it "breathes" versus time, then the cumulative stats. will be smeared. Instantaneously, though, the distribution will be bi-modal.

@ Simon,

Well, there you go. If you can detect a signature of the bid-ask spread in the trade data, then perhaps it's enough information to help in a trading strategy. And if you can't detect it, then the strategy could just stay out of securities with unknown estimated bid-ask until the spread can be detected above a level of statistical significance. Another thought is that over short enough time scales, perhaps the drift in the mid-point (average of bid-ask prices) could be modeled as a low-order polynomial (linear or possibly quadratic), over N minutes intra-day.

Regarding thinly traded securities, I don't understand why, at an instant of time, there would always be both bid and ask prices. Say I want to sell 100 shares of PIG (Pork Industries Gigante), my massive Midwestern farming enterprise. If nobody is interested in PIG, then when I put the lot out for sale, there will be only an ask price, right? Or does a bid price always magically appear since there will always be some algorithm following PIG that will automatically pick a bid of \$0.

Generally, because designated market makers are required to make two-sided markets. I don't know if there are any stocks without designated market makers or specialists, it's been a while since I was involved in that stuff, and I think some rules changed after the flash crash.

Just a FYI: I tried something similar some time ago, see here: https://github.com/harpone/PyArb

The holding times I used were very short and profits were pretty much killed by transaction costs... I think it's quite a bit simpler to implement than the other model in regards to estimation etc. Also IMO the model is way too simple, and in fact it might be smarter to model the returns instead of (log) prices (there are additional subtleties in estimating the price time series process).

Never seemed to have enough time to try that out in Quantopian, but if anyone gets something out of it, let me know :)

Hi Pravin,
I was wondering if you tried to implement this nice algo with CVXPY to scale with bigger stock universe ?

Hi Mat,

I wasted several months on this algorithm only to realize that it relies on covariance matrix in level data (non-stationary). I have since then discarded this approach to look for alternative techniques. Will port it to CVXPY if time permits.

Best regards,
Pravin

@Mat, if you read page 3 of this paper, you will see that it assumes that asset prices follow a stationary VAR process which is not at all reflective of true nature of stock prices. It could work on something like implied volatility but Quantopian does not have option data.

Thanks for your answer Pravin, maybe we can add some regime switching model to this algorithm like HMM or other based on volatility to detect changes in series behavior ?
What is your last research for finding stationary asset prices ?
Regards,

@Matt, If could find a way to compute stationary covariance matrices from non-stationary prices, it will be revolutionary. My current research is much simpler based on return series which are at least weakly stationary.


text:
[This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

### TL;DR

• C++ templates and function overloading are incompatible with R’s C API, so polymorphism must be achieved via run-time dispatch, handled explicitly by the programmer.
• The traditional technique for operating on SEXP objects in a generic manner entails a great deal of boilerplate code, which can be unsightly, unmaintainable, and error-prone.
• The desire to provide polymorphic functions which operate on vectors and matrices is common enough that Rcpp provides the utility macros RCPP_RETURN_VECTOR and RCPP_RETURN_MATRIX to simplify the process.
• Subsequently, these macros were extended to handle an (essentially) arbitrary number of arguments, provided that a C++11 compiler is used.

### Background

To motivate a discussion of polymorphic functions, imagine that we desire a function (ends) which, given an input vector x and an integer n, returns a vector containing the first and last n elements of x concatenated. Furthermore, we require ends to be a single interface which is capable of handling multiple types of input vectors (integers, floating point values, strings, etc.), rather than having a separate function for each case. How can this be achieved?

#### R Implementation

A naïve implementation in R might look something like this:

[1] 1 2 3 4 6 7 8 9

[1] "a" "b" "c" "x" "y" "z"

[1] -0.560476 -0.230177  0.701356 -0.472791

The simple function above demonstates a key feature of many dynamically-typed programming languages, one which has undoubtably been a significant factor in their rise to popularity: the ability to write generic code with little-to-no additional effort on the part of the developer. Without getting into a discussion of the pros and cons of static vs. dynamic typing, it is evident that being able to dispatch a single function generically on multiple object types, as opposed to, e.g. having to manage separate impementations of ends for each vector type, helps us to write more concise, expressive code. Being an article about Rcpp, however, the story does not end here, and we consider how this problem might be approached in C++, which has a much more strict type system than R.

#### C++ Implementation(s)

For simplicity, we begin by considering solutions in the context of a “pure” (re: not called from R) C++ program. Eschewing more complicated tactics involving run-time dispatch (virtual functions, etc.), the C++ language provides us with two straightforward methods of achieving this at compile time:

1. Function Overloading (ad hoc polymorphism)
2. Templates (parametric polymorphism)

The first case can be demonstrated as follows:

Although the above program meets our criteria, the code duplication is profound. Being seasoned C++ programmers, we recognize this as a textbook use case for templates and refactor accordingly:

This approach is much more maintainable as we have a single implementation of ends rather than one implementation per typedef. With this in hand, we now look to make our C++ version of ends callable from R via Rcpp.

#### Rcpp Implementation (First Attempt)

Many people, myself included, have attempted some variation of the following at one point or another:

Sadly this does not work: magical as Rcpp attributes may be, there are limits to what they can do, and at least for the time being, translating C++ template functions into something compatible with R’s C API is out of the question. Similarly, the first C++ approach from earlier is also not viable, as the C programming language does not support function overloading. In fact, C does not support any flavor of type-safe static polymorphism, meaning that our generic function must be implemented through run-time polymorphism, as touched on in Kevin Ushey’s Gallery article Dynamic Wrapping and Recursion with Rcpp.

#### Rcpp Implementation (Second Attempt)

Armed with the almighty TYPEOF macro and a SEXPTYPE cheatsheat, we modify the template code like so:

[1] 1 2 3 4 6 7 8 9

[1] "a" "b" "c" "x" "y" "z"

[1] -1.067824 -0.217975 -0.305963 -0.380471

Warning in ends(list()): Invalid SEXPTYPE 19 (VECSXP).

NULL

Some key remarks:

1. Following the ubiquitous Rcpp idiom, we have converted our ends template to use an integer parameter instead of a type parameter. This is a crucial point, and later on, we will exploit it to our benefit.
2. The template implementation is wrapped in a namespace in order to avoid a naming conflict; this is a personal preference but not strictly necessary. Alternatively, we could get rid of the namespace and rename either the template function or the exported function (or both).
3. We use the opaque type SEXP for our input / output vector since we need a single input / output type. In this particular situation, replacing SEXP with the Rcpp type RObject would also be suitable as it is a generic class capable of representing any SEXP type.
4. Since we have used an opaque type for our input vector, we must cast it to the appropriate Rcpp::Vector type accordingly within each case label. (For further reference, the list of vector aliases can be found here). Finally, we could dress each return value in Rcpp::wrap to convert the Rcpp::Vector to a SEXP, but it isn’t necessary because Rcpp attributes will do this automatically (if possible).

At this point we have a polymorphic function, written in C++, and callable from R. But that switch statement sure is an eyesore, and it will need to be implemented every time we wish to export a generic function to R. Aesthetics aside, a more pressing concern is that boilerplate such as this increases the likelihood of introducing bugs into our codebase – and since we are leveraging run-time dispatch, these bugs will not be caught by the compiler. For example, there is nothing to prevent this from compiling:

// ...
case INTSXP: {
return impl::ends(as<CharacterVector>(x), n);
}
// ...

In our particular case, such mistakes likely would not be too disastrous, but it should not be difficult to see how situations like this can put you (or a user of your library!) on the fast track to segfault.

### Obligatory Remark on Macro Safety

The C preprocessor is undeniably one of the more controversial aspects of the C++ programming language, as its utility as a metaprogramming tool is rivaled only by its potential for abuse. A proper discussion of the various pitfalls associated with C-style macros is well beyond the scope of this article, so the reader is encouraged explore this topic on their own. On the bright side, the particular macros that we will be discussing are sufficiently complex and limited in scope that misuse is much more likely to result in a compiler error than a silent bug, so practically speaking, one can expect a fair bit of return for relatively little risk.

### Synopsis

At a high level, we summarize the RCPP_RETURN macros as follows:

• There are two separate macros for dealing with vectors and matrices, RCPP_RETURN_VECTOR and RCPP_RETURN_MATRIX, respectively.
• In either case, code is generated for the following SEXPTYPEs:
• INTSXP (integers)
• REALSXP (numerics)
• RAWSXP (raw bits)
• LGLSXP (logicals)
• CPLXSXP (complex numbers)
• STRSXP (characters / strings)
• VECSXP (lists)
• EXPRSXP (expressions)
• In C++98 mode, each macro accepts two arguments:
1. A template function
2. A SEXP object
• In C++11 mode (or higher), each macro additionally accepts zero or more arguments which are forwarded to the template function.

Finally, the template function must meet the following criteria:

• It is templated on a single, integer parameter.
• In the C++98 case, it accepts a single SEXP (or something convertible to SEXP) argument.
• In the C++11 case, it may accept more than one argument, but the first argument is subject to the previous constraint.

Examining our templated impl::ends function from the previous section, we see that it meets the first requirement, but fails the second, due to its second parameter n. Before exploring how ends might be adapted to meet the (C++98) template requirements, it will be helpful demonstrate correct usage with a few simple examples.

### Fixed Return Type

We consider two situations where our input type is generic, but our output type is fixed:

1. Determining the length (number of elements) of a vector, in which an int is always returned.
2. Determining the dimensions (number of rows and number of columns) of a matrix, in which a length-two IntegerVector is always returned.

First, our len function:

(Note that we omit the return keyword, as it is part of the macro definition.) Testing this out on the various supported vector types:

[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Similarly, creating a generic function that determines the dimensions of an input matrix is trivial:

And checking this against base::dim,

[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

everything seems to be in order.

It’s worth pointing out that, for various reasons, it is possible to pass a matrix object to an Rcpp function which calls RCPP_RETURN_VECTOR:

[1] 9

[1] 9

Although this is sensible in the case of len – and even saves us from implementing a matrix-specific version – there may be situations where this behavior is undesirable. To distinguish between the two object types we can rely on the API function Rf_isMatrix:

[1] 9

<Rcpp::exception in len2(matrix(1:9, 3)): matrix objects not supported.>

We don’t have to worry about the opposite scenario, as this is already handled within Rcpp library code:

<Rcpp::not_a_matrix in dims(1:5): Not a matrix.>

### Generic Return Type

In many cases our return type will correspond to our input type. For example, exposing the Rcpp sugar function rev is trivial:

[1] 5 4 3 2 1

[[1]]
[1] 5+2i

[[2]]
[1] 4+2i

[[3]]
[1] 3+2i

[[4]]
[1] 2+2i

[[5]]
[1] 1+2i

[1] "edcba"

As a slightly more complex example, suppose we would like to write a function to sort matrices which preserves the dimensions of the input, since base::sort falls short of the latter stipulation:

[1] 1 2 3 4 5 6 7 8 9

There are two obstacles we need to overcome:

1. The Matrix class does not implement its own sort method. However, since Matrix inherits from Vector, we can sort the matrix as a Vector and construct the result from this sorted data with the appropriate dimensions.
2. As noted previously, the RCPP_RETURN macros will generate code to handle exactly 8 SEXPTYPEs; no less, no more. Some functions, like Vector::sort, are not implemented for all eight of these types, so in order to avoid a compilation error, we need to add template specializations.

With this in mind, we have the following implementation of msort:

Note that elements will be sorted in column-major order since we filled our result using this constructor. We can verify that msort works as intended by checking a few test cases:

     [,1] [,2] [,3]
[1,]    1    7    4
[2,]    3    9    6
[3,]    5    2    8

     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

[1] 1 2 3 4 5 6 7 8 9

     [,1] [,2]
[1,] "a"  "y"
[2,] "c"  "b"
[3,] "z"  "x"

     [,1] [,2]
[1,] "a"  "x"
[2,] "b"  "y"
[3,] "c"  "z"

[1] "a" "b" "c" "x" "y" "z"

List of 9
$: int 1$ : int 2
$: int 3$ : int 4
$: int 5$ : int 6
$: int 7$ : int 8
\$ : int 9
- attr(*, "dim")= int [1:2] 3 3

<Rcpp::exception in msort(x): sort not allowed for lists.>

<simpleError in sort.int(x, na.last = na.last, decreasing = decreasing, ...): 'x' must be atomic>

### Revisiting the ‘ends’ Function

Having familiarized ourselves with basic usage of the RCPP_RETURN macros, we can return to the problem of implementing our ends function with RCPP_RETURN_VECTOR. Just to recap the situation, the template function passed to the macro must meet the following two criteria in C++98 mode:

1. It is templated on a single, integer parameter (representing the Vector type).
2. It accepts a single SEXP (or convertible to SEXP) argument.

Currently ends has the signature

meaning that the first criterion is met, but the second is not. In order preserve the functionality provided by the int parameter, we effectively need to generate a new template function which has access to the user-provided value at run-time, but without passing it as a function parameter.

The technique we are looking for is called partial function application, and it can be implemented using one of my favorite C++ tools: the functor. Contrary to typical functor usage, however, our implementation features a slight twist: rather than using a template class with a non-template function call operator, as is the case with std::greater, etc., we are going to make operator() a template itself:

Not bad, right? All in all, the changes are fairly minor:

• The function body of Ends::operator() is identical to that of impl::ends.
• n is now a private data member rather than a function parameter, which gets initialized in the constructor.
• Instead of passing a free-standing template function to RCPP_RETURN_VECTOR, we pass the expression Ends(n), where n is supplied at run-time from the R session. In turn, the macro will invoke Ends::operator() on the SEXP (RObject, in our case), using the specified n value.

We can demonstrate this on various test cases:

[1] 1 2 3 4 6 7 8 9

[1] "a" "b" "c" "x" "y" "z"

[1] -0.694707 -0.207917  0.123854  0.215942

### A Modern Alternative

As alluded to earlier, a more modern compiler (supporting C++11 or later) will free us from the “single SEXP argument” restriction, which means that we no longer have to move additional parameters into a function object. Here is ends re-implemented using the C++11 version of RCPP_RETURN_VECTOR (note the // [[Rcpp::plugins(cpp11)]] attribute declaration):

[1] 1 2 3 4 6 7 8 9

[1] "a" "b" "c" "x" "y" "z"

[1]  0.379639 -0.502323  0.181303 -0.138891

The current definition of RCPP_RETURN_VECTOR and RCPP_RETURN_MATRIX allows for up to 24 arguments to be passed; although in principal, the true upper bound depends on your compiler’s implementation of the __VA_ARGS__ macro, which is likely greater than 24. Having said this, if you find yourself trying to pass around more than 3 or 4 parameters at once, it’s probably time to do some refactoring.

To leave a comment for the author, please follow the link and comment on their blog: Rcpp Gallery.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

# Never miss an update! Subscribe to R-bloggers to receive e-mails with the latest R posts.(You will not see this message again.)

Click here to close (This popup will not appear again)


text:
## Asterisms in LaTeX

February 22, 2011

Another $\LaTeX$ trick this week: asterisms. Asterisms are typographic devices used to separate or to call attention to a piece of text without resorting to using a section or chapter. While its usage in modern English (and French, for that matter) is now rare, it is still an interesting device to structure text.

The usual manifestation of an asterism is the ⁂ symbol (U+2042) but the single character asterism may not be suitable for all occasion. In HTML, and in these posts, you have seen it quite often. With $\LaTeX$, getting an decent-looking asterism is not that easy.

## Epigraphs in LaTeX

January 18, 2011

There are times when part of the message, the gist, must be communicated to the reader in an out-of-band fashion, so to speak. One way of doing this is to use an epigraph to open a chapter or section, carefully chosen to convey the intended message but in the voice of another author (self-epigraphs are of very bad taste in my opinion).

$\LaTeX$ is the preferred document preparation system of computer scientists, physicists, and mathematicians and if you intend to follow a career into the academia, it’s pretty much unavoidable. One day, you’ll have to learn $\LaTeX$. The thing is, $\LaTeX$ is pretty much like C++: it can do just about anything, but it’s not going to help you do it. You have to rely on the innumerable packages or, if you really can’t find what you need, you can code it yourself. Let us have a look on how to code an epigraph macro in $\LaTeX$.

## Suggested Reading:The LATEX Companion

September 4, 2009

Frank Mittelbach, Michel Goossens, Johannes Braams, David Carlisle, Chris Rowley — The LATEX Companion — 2nd ed, Addison Wesley, 2006, 1090 pp. ISBN 0-201-36299-6


text:
# Find the fastest path through this Brazilian concretist poem

One day, a small jet was flying over the concretist poem "Velocidade" (Speed), by the Brazilian poet Ronaldo Azeredo.

In order to cover every vowel present in this poem, it goes in straight line from one letter to another either orthogonally or diagonally adjacent, never making a turn sharper than 45 degrees*.

Assuming the letters are like points in a square lattice, and the jet travels at constant speed, what is the fastest way this can be achieved by the pilot?

Notice that no restrictions were placed in where the jet should start, except that its position goes from one letter to another. Also, the jet is not allowed to leave the matrix of letters.

* Measured by the external angle, i.e. deviation from the original path, so that a 0 degree turn leads to the same direction, while a 180 degree turn makes the jet go in the reverse.

• ...Where does the pilot begin? – warspyking Sep 21 '15 at 13:26
• Can the plane leave the square once it has entered? – LogicianWithAHat Sep 21 '15 at 14:34
• You'd have to end on the E at the bottom-right corner, because otherwise there's no way for you to turn out of it. – Joe Z. Sep 21 '15 at 15:37
• OP doesn't say it's not allowed to use vowels several times. (This would not be solvable otherwise.) – Tom K. Sep 21 '15 at 16:07
• Do you mean "sharper than 135 degrees"? All your answers seem to heed this constraint. – djechlin Sep 21 '15 at 18:31

Here's a solution with

54 steps.

$38 \space horizontal/vertical + 16 \space diagonal => length \space 60.627$
Thanks to f'' for length optimization.

A brute-forced solution with

47 steps.

$28 \space horizontal/vertical + 19 \space diagonal => length \space 54.87$

• With three extra vowel touches in this one, do you think hypothetically there should be a solution with 51, or is that not guaranteed? – lorimer Sep 21 '15 at 18:13
• @lorimer I think the extra vowels are no more different than the consonants. It's possible there is a better solution. – Sleafar Sep 21 '15 at 18:18
• It looks like you have a 90 degree turn at the bottommost 'O' – Sconibulus Sep 21 '15 at 20:46
• You can save 4(sqrt2-1) by changing the triple crossing to three vertical lines. – f'' Sep 21 '15 at 20:57
• You can also save another 2(sqrt2-1) by changing the LC/OO crossing in the lower middle to two horizontal segments. – f'' Sep 21 '15 at 21:00

OK. Touching every vowel once, using coordinates A-J horizontal and 0-9 vertically, with the origin in the top left:

J9->J1, I0, H0->A7, A8, B9->H9
I9->I2, H1, G1->B6, B7, C8->G8
H7->H3, G2, F2->C5, C6, D7->F7
G6->G4, F3->D3, C4, C5, D6->E6, F5

That's 71 steps. It's probably not optimized perfectly, so I'll keep looking, but it's at least a starting point.

EDIT: This one has 59 but touches some vowels more than once. Counting as distance (cf. Sleafar's comment) it would be 68.11 units with 22 diagonal moves.

J9->J1, I0, H0, G1, G2, H3, I4->I8
H9->B9, A8, A7, B6, C6->E8->G8, H7
H7->H5, G4, F4->C7, C8, D9, E9, F8
F7->I4->I2, H1, G1, F2->F5, E6

This one doesn't waste as much time by avoiding U-turns, but instead grabs what it can that's nearby.

ETA3: nevermind, the edit would have missed one 'E'. i'll quit second-guessing myself now

I think I am missing something obvious, since the solution couldn't be that easy

16 diagonal + 16 vertical/horizontal steps.

• "never making a turn sharper than 45 degrees" – Tryth Sep 22 '15 at 3:11
• where? I am brain-dead here, and can see only 90 degrees and 135 degrees – kushj Sep 22 '15 at 3:25
• Since this point has been brought up before, I edited my post to clarify it – MathET Sep 22 '15 at 3:41


text:
# Comparing probabilities of drawing balls of certain color, with and without replacement

A bag contains 4 white 2 black 3 red balls. A ball is drawn from the bag 1-by-1.
(i) Find the probability that $5^{th}$ ball is red. Consider both with and without replacement.
(ii) 3^rd and 7^th are red. Find the probability that $9^{th}$ ball is black. Consider both with and without replacement.

### (i)

(a) with replacement. $$P=\frac39=\frac13$$

(b) without replacement. Let's make an array of 9 elements where i^th element denote which ball is taken out in i^th turn, then we can have red in 5^th position with probability:$$P=\frac{\binom31\binom82\binom622!2!4!}{9!}=\frac13$$ [select one red ball to put in 5^th position, then select positions on array for black and remaining red then arrange them divided by total arrangements]

### (ii)

(a) with replacement. $$P=\frac29$$

(b) without replacement. Let's make an array of 9 elements where i^th element denote which ball is taken out in i^th turn, then we have red in 3rd and 7th position red,we can have black in 9^th position with probability:$$P=\frac{\binom32\binom212!6!2!}{\binom322!7!}=\frac27$$ [select two red balls for 3 and 7 position then select one black ball for last position and then select position for red and remaining black balls on array then arrange]

Why did drawing with or without replacement make no difference in the first and second can be obtained by considering 3^rd and 7^th balls were never there?

• Not sure why you said we must read your solution so very carefully. But do see if my answer addresses your real question. =) – user21820 Nov 29 '14 at 5:54
• @user21820 it would take me some time, thanks – RE60K Nov 29 '14 at 5:56
• As turkeyhundt pointed out, the probability for (ii) with replacement is 2/9, not 2/7. It is 2/7 without replacement. – user21820 Nov 29 '14 at 6:54
• @user21820 edited. – RE60K Nov 29 '14 at 8:14

We can consider each configuration for without replacement as a sequence of balls. Then permuting the sequence in any way gives another configuration with equal probability because every sequence has equal probability. Each permutation can be used to give a bijection (one-to-one correspondence) between a set of configurations with another that differ simply by that permutation.

For example, in the first question you want #(sequences with 5th ball red)/#(sequences), and #(sequences with 5th ball red) = #(sequences with 1st ball red) by the bijection that switches the 1st and 5th ball in a sequence. Check that every sequence with 5th ball red corresponds to exactly one sequence with 1st ball red, and vice versa. This gives the answer since #(sequences with 1st ball red)/#(sequences) = 3/9 because there are no restrictions on the later balls.

For the second question, to do it rigorously we choose some permutation that sends the 3rd ball to position 1, the 7th ball to position 2 and the 9th ball to position 3. Then we get #(sequences with 3rd and 7th ball red and 9th ball black) = #(sequences with 1st and 2nd ball red and 3rd ball black), and #(sequences with 3rd and 7th ball red) = #(sequences with 1st and 2nd ball red). So the answer is #(sequences with 1st and 2nd ball red and 3rd ball black)/#(sequences with 1st and 2nd ball red) = 2/7 since the first two balls already use up two red and after drawing the first two we are left with 7 balls of which 2 are black.

• i can't understand "If we switch the 1st and k-th ball in the sequence, it is in both your questions a bijection between valid configurations." – RE60K Nov 29 '14 at 6:08
• I should have been more precise. In your first question, you are interested in the 5th ball. Given any sequence drawn, we can switch the 1st and 5th ball. That results in another possible drawn sequence, and this switching process is a bijection between the set of drawn sequences with 5th ball of interest and the set of drawn sequences with 1st ball of interest. Similarly for the second question, but we have to exclude the used balls. I can expand my answer if it's not clear. – user21820 Nov 29 '14 at 6:21
• you are very close to explaining me i understand exchnage but what do you mean by bijection, bijective functions which are invertible? – RE60K Nov 29 '14 at 6:27
• A bijection is a exact correspondence between two sets, meaning that each element of one set corresponds to exactly one element of the other set. I'll expand my answer to make it more rigorous since it's not so trivial for your second question. – user21820 Nov 29 '14 at 6:32
• @eaxdpiotnyeantial: See if the above makes sense now. It should be as rigorous as possible without clutter. Note the important points are that every sequence is equally likely, which you need to prove separately as a lemma, and that a permutation is really a bijection as claimed, which is intuitive but should also be given a proof. – user21820 Nov 29 '14 at 6:45

Without replacement is a word. So consider strings of the form $\{R, B, G\}^{*}$. That is, strings consisting only of letters $R$, $B$, and $G$. If we do replace, there are $3^{n}$ strings of length $n$ over this alphabet. That is, we have an inexhaustable supply of each character.

If we do not replace, we have $3!$ strings of length $3$ (permute the letters), $3!$ strings of length $2$, and $3$ strings of length $1$. We cannot form a string of length $4$.

So consider:

A bag contains 4 white 2 black 3 red balls. A ball is drawn from the bag 1-by-1. (i) Find the probability that 5th ball is red.Consider both w and w/o replacement.

If we do replace, the odds of getting a red ball on the $nth$ character are the same as for the first character. So the red ball on $5$ is $\frac{3}{9}$.

Combinatorially, if we have a string of length $n \geq 5$, we know a character will be drawn for each slot. Let's label each ball to make the visualization a little easier, so we have $9$ characters. We are interested in the 5th being one of $r_{1}, r_{2}, r_{3}$. So there are $9^{n-1} * 3$ such strings with red on character $5$, and $9^{n}$ such strings of length $n$. Do you see how I got that combinatorially based on my explanation before your problem?

Now if we do not replace, the explanation you described in your solution is spot on. This is a permutation problem, not a word problem.

• what do you mean by "word problem"? – RE60K Nov 29 '14 at 6:41
• A "word" problem is where we have an alphabet like $\{R, G, B\}$ and we construct strings or "words" over that alphabet. – ml0105 Nov 29 '14 at 6:43

Consider this simpler problem, conceptually, as an example.

What if you have a blue, green, red, yellow, purple and black marble. What if you draw 3 marbles without replacement. What is the probability the 3rd one is black? Well, no marble has an inherent advantage over the others, so intuitively, it must be a one sixth chance.

Now, do the same problem, but WITH replacement. Again, no marble has an inherent advantage or disadvantage to being chosen, so there must be a $\frac{1}{6}$ chance of the 3rd one being black.

• this explains the first part very well, could you expand your answer further to include second part as well? – RE60K Nov 29 '14 at 6:40
• Hmmm, let me take a look. I'm not convinced WITH REPLACEMENT, that the probability would be $\frac{2}{7}$. It seems like it would be $\frac{2}{9}$ if you are replacing your ball each time. Am I reading the problem correctly? – turkeyhundt Nov 29 '14 at 6:46
• Also, one thing that might help on the WITHOUT REPLACEMENT is, instead of imagining them jumbled up in a bag, imagine them hidden, but queued up in a pipe or half pipe and they are taken from the end one at a time. Once they are shuffled, the outcome is set. So regardless of what is chosen ahead of it, each ball's fate is already sealed. – turkeyhundt Nov 29 '14 at 6:51
• You're right. I didn't notice that with replacement the known ball colours become irrelevant. Then I guess the real question is why we can sort of see what happens if we have taken out the known balls. Your idea is indeed leads to the method I outline in my answer. – user21820 Nov 29 '14 at 6:52


text:
# nLab Einstein-Hilbert action

### Context

#### Gravity

gravity, supergravity

# Contents

## Idea

The Einstein-Hilbert action is the action functional that defines the dynamics of gravity in general relativity. It is a canonical invariant of pseudo-Riemannian manifolds.

## Definition

For $\left(X,g\right)$ a Riemannian manifold or pseudo-Riemannian manifold, its vacuum Einstein-Hilbert action is the number

$\left(X,g\right)↦{\int }_{X}R\left(g\right)\mathrm{dvol}\left(g\right)\phantom{\rule{thinmathspace}{0ex}},$(X,g) \mapsto \int_X R(g) dvol(g) \,,

where

• $R\left(g\right):X\to ℝ$ is the scalar curvature function of the metric $g$;

• $\mathrm{dvol}\left(g\right)$ is the volume form defined by the metric

• ${\int }_{X}$ is integration of differential forms over $X$.

If the metric is instead encoded in terms of an Poincare group-connection $\nabla$ (the first-order formulation of gravity) then (for $\mathrm{dim}X=4$ and assuming for simplicity that the underlying bundle is trivial) this is equivalently

$\left(X,\nabla \right)↦{\int }_{X}⟨R\wedge e\wedge e⟩\phantom{\rule{thinmathspace}{0ex}},$(X,\nabla) \mapsto \int_X \langle R \wedge e \wedge e \rangle \,,

where now

• $R$ denotes the curvature 2-form of the orthogonal group-connection part of $\nabla$;

• $e$ denotes the vielbein, (the translation part of the connection 1-form);

• $⟨-⟩$ is the invariant polynomial which in terms of the canonical basis has components $\left({ϵ}_{abcd}\right)$.

This extends to an action functional for gravity coupled to “matter” (which in the context of general relativity conventionally means every other field, for instance the electromagnetic field) by simply adding the matter action on a given gravitational background

$\left(X,g,\varphi \right)↦{\int }_{X}\left({L}_{\mathrm{EH}}\left(g\right)+{L}_{\varphi }\left(g\right)\right)d\mathrm{vol}\left(g\right)\phantom{\rule{thinmathspace}{0ex}}.$(X,g,\phi) \mapsto \int_X (L_{EH}(g) + L_\phi(g)) d vol(g) \,.

The critical points of the Einstein-Hilbert action define the physically realized spacetimes. This are Einstein's equations (the Euler-Lagrange equations of this action)

$\delta {S}_{\mathrm{EH}}\left(g\right)=0\phantom{\rule{thinmathspace}{0ex}}.$\delta S_{EH}(g) = 0 \,.

(…)

standard model of particle physics and cosmology

theory:Einstein-Yang-Mills-Dirac-Higgs
gravityelectroweak and strong nuclear forcefermionic matterscalar field
field content:vielbein field $e$principal connection $\nabla$spinor $\psi$scalar field $H$
Lagrangian:scalar curvature densityfield strength squaredDirac operator component densityfield strength squared + potential density
$L=$$R\left(e\right)\mathrm{vol}\left(e\right)+$$⟨{F}_{\nabla }\wedge {\star }_{e}{F}_{\nabla }⟩+$$\left(\psi ,{D}_{\left(e,\nabla \right)}\psi \right)\mathrm{vol}\left(e\right)+$$\nabla \overline{H}\wedge {\star }_{e}\nabla H+\left(\lambda {\mid H\mid }^{4}-{\mu }^{2}{\mid H\mid }^{2}\right)\mathrm{vol}\left(e\right)$

## References

Section Prequantum gauge theory and Gravity at

Revised on January 14, 2013 18:20:46 by Urs Schreiber (203.116.137.162)


text:
# Can a primitive root-permutation of $A=\{1, 2, \ldots, p-1 \}$ be a cycle of length $p-1$ only for finitely many $p$?

Let $p$ be an odd prime, $g$ a primitive root of $p$ and $A=\{1, 2, \ldots, p-1 \}$.

Obviously, $\sigma_g(p)=\begin{pmatrix} 1 & 2 & \ldots & {p-1} \\ g^1\pmod{p} & g^2\pmod{p} & \ldots & g^{p-1}\pmod{p} \end{pmatrix}$ is a permutation of $A$.

I observed that "almost always" $\sigma$ is a product of cycles whose length is strictly less than $p-1$.
If $p=3$ and $g=2$, then $\sigma_2(3)=\begin{pmatrix} 1 & 2 \\ 2 & 1 \end{pmatrix}=(1\quad2)$ is a cycle of length $3-1=2$.

If $p=5$ and $g=3$, then $\sigma_3(5)=\begin{pmatrix} 1 & 2 & 3 & 4 \\ 3 & 4 & 2 & 1 \end{pmatrix}=(1\quad 3\quad 2\quad 4)$ is a cycle of length $5-1=4$.
But these two cases seem to be just exceptions.
For other (small) values of $p$ and every primitive root $g$ of $p$, there are no cycles of length $p-1$.

I can prove that $g=2$ and $g=\frac{p+1}{2}$ cannot produce such a big cycle, but I do not have any idea how to attack the problem for other values of $g$.
I conjecture that for every $g$ there is not such a big cycle (for sufficiently large $p$)

My question is:

Do only finitely many primes $p$ exist, such that for some primitive root $g$ the permutation $\sigma_g(p)=\begin{pmatrix} 1 & 2 & \ldots & {p-1} \\ g^1\pmod{p} & g^2\pmod{p} & \ldots & g^{p-1}\pmod{p} \end{pmatrix}$ is a cycle?

• Apparently, there are a lot of counterexamples. Smallest ones are $p = 23, g = 20$; $p = 41, g = 6$; $p = 59, g = 39$; $p = 61, g = 10$; $p = 107, g = 94$. They don't stop appearing as $p$ grows either (and there are $p$-s, for which more than one $g$ is suitable). So it seems that the answer is more probably "No" than "Yes". – Kaban-5 Jul 29 '18 at 22:04
• @Kaban-5 what are the first few $p$ for which there is more than one $g$ leading to a $(p-1)$-cycle, and what are all such $g$ for those $p$? – KConrad Jul 30 '18 at 11:42
• Two possible $g$'s: $p = 587, g = 150, 375$; $p = 751, g = 240, 263$; $p = 809, g = 265, 750$; $p = 811, g = 113, 165$. Three: $p = 1889, g = 479, 859, 1248$; $p = 2267, g = 61, 976, 2010$; $p = 2699, g = 427, 639, 1256$; $p = 3491, g = 1472, 1626, 1833$. Four: $p = 5417, g = 1319, 1566, 3486, 5290$; $p = 7691, g = 558, 1434, 6244, 6760$. Five and more does not happen for $p \leqslant 10^4$. I believe that this should happen eventually anyway because heuristics suggest that, at least assuming that the number of Sophie Germain primes is infinite (probably unnecessary). – Kaban-5 Jul 30 '18 at 16:21

As pointed out in the comments, primitive roots with a single cycle appear to be rather common. The standard heuristic argument suggests that there are infinitely many such $p$, and a bit more.
The share of $N$-cycles among all permutations on $N$ symbols is $1/N$. By the Borel-Cantelli lemma, given an infinite sequence $\{\sigma_i\}$ of independent uniformly distributed random permutations on $N_i$ symbols, if $\sum 1/N_i$ diverges then with probability one infinitely many $\sigma_i$ constitute a single cycle. Since the series of reciprocal primes diverges, we may expect that given a random sequence $\{g_i\}$ of primitive roots mod $p_i$ over all primes, with probability one there are infinitely many $\sigma_i=\sigma_{g_i}(p_i)$ that constitute a single cycle. Of course, since uniform distribution and independence are only heuristic, this is not a proof!


text:
## Friday, May 27, 2011

### Ellipse detection in image by using Hough transform

How we can detect ellipses in images ? One way is to use Hough transform. I will use Hough transform algorithm variant created by Yonghong Xie and Qiang Ji. That algorithm pseudo-code:

1. Store all edge pixels in a one dimensional array.
2. Clear the accumulator array.
3. For each pixel (x1, y1), carry out the following steps from (4) to (14).
4. For each other pixel (x2, y2), if the distance between
(x1, y1) and (x2, y2) is greater than the required least
distance for a pair of pixels to be considered then
carry out the following steps from (5) to (14).
5. From the pair of pixels (x1, y1) and (x2, y2) calculate the center,
orientation and length of major axis for the assumed ellipse.
6. For each third pixel (x, y), if the distance between
(x, y) and (x0, y0) is greater than the required least
distance for a pair of pixels to be considered then
carry out the following steps from (7) to (9).
7. Calculate the length of minor axis.
8. Increment the accumulator for this length of minor axis by 1.
9. Loop until all pixels are computed for this pair of pixels.
10. Find the maximum element in accumulator array. The
related length is the possible length of minor axis
for assumed ellipse. If the vote is greater than the
required least number for assumed ellipse, one
ellipse is detected.
11. Output ellipse parameters.
12. Remove the pixels on the detected ellipse from edge pixel array.
13. Clear accumulator array.
14. Loop until all pairs of pixels are computed.

Proof-of-concept algorithm implementation in Python:
import sys
from PIL import Image,ImageFilter, ImageDraw
from math import *

# some global constants
EL_COVERAGE_RATIO = 0.9
EL_VERIFICATION_DISTANCE = 1.
EL_PATH_POINTS = 51
MIN_MINOR_FREQUENCY = 30
MIN_HALF_MAJOR = 8
MIN_HALF_MINOR = 6

def distance(p1,p2):
x1,y1 = p1
x2,y2 = p2
return sqrt((x1-x2)**2 + (y1-y2)**2)

def nonnegative(v):
return v if v >= 0 else 0

def parametricEllipse(center, a, b, angle):
xc,yc = center
elx = lambda t: xc + a * cos(t) * cos(angle) - b * sin(t) * sin(angle)
ely = lambda t: yc + a * cos(t) * sin(angle) + b * sin(t) * cos(angle)
return [(int(elx(2.*pi*x/float(EL_PATH_POINTS-1))),int(ely(2.*pi*x/float(EL_PATH_POINTS-1)))) for x in range(EL_PATH_POINTS)]

assert len(sys.argv)!=2, "missing input and/or output file !"

im = Image.open(sys.argv[1])
width, height = im.size
io = Image.new('RGB',(width, height),(255,255,255))
draw = ImageDraw.Draw(io)

# converting image to grayscale
im = im.convert('L')

# detecting edge pixels
im = im.filter(ImageFilter.FIND_EDGES)

# converting to binary image
im = im.convert('1')
pxy = []

# extracting binary pixels coordinates
for x in range(width):
for y in range(height):
if pixels[x,y]==255:
pxy.append((x,y))

# applying Hough transform for ellipses detection.
# algorithm is taken from this paper:
cIx = -1
colors = [(255,0,0),(0,200,0),(0,0,255)]
for x1,y1 in pxy:
for x2,y2 in pxy:
bbins = {}
dist = distance((x1,y1),(x2,y2))
if dist >= 2*MIN_HALF_MAJOR:
cent = ((x1+x2)/2.,(y1+y2)/2.)
a = dist/2. # semi-length of major axis
alfa = atan2((y2 - y1),(x2 - x1))
for rx,ry in pxy:
d = distance((rx,ry),cent)
if d >= MIN_HALF_MINOR:
f = distance((rx,ry),(x2,y2))
cost = (a**2. + d**2. - f**2.)/(0.00001+2.*a*d)
b = sqrt(nonnegative((a**2. * d**2. * (1.-cost**2.))/(0.00001 + a**2. - d**2. * cost**2.)))  # semi-length of minor axis
b = int(b)
if bbins.has_key(b):
bbins[b]+=1
elif b > 0:
bbins[b]=1
bbins_rev = dict([(v,k) for k,v in bbins.iteritems()])
max_freq = max(bbins_rev.keys())
bmax = bbins_rev[max_freq]
# Did we found probable ellipse ?
if max_freq >= MIN_MINOR_FREQUENCY and alfa >=0.0 and bmax >= MIN_HALF_MINOR:
elData = parametricEllipse(cent, a, bmax, alfa)
supported = []
supportRatio = 0.0
# counting how much pixels lies on ellipse path
for i in range(EL_PATH_POINTS):
elx,ely = elData[i]
for x,y in pxy:
if distance((elx,ely),(x,y)) <= EL_VERIFICATION_DISTANCE:
supported.append((x,y))
supportRatio += 1./float(EL_PATH_POINTS)
supported = list(set(supported))
# if number of pixels on ellipse path is big enough
if supportRatio >= EL_COVERAGE_RATIO:
cIx = (cIx+1)%3
print "coverage %.2f" % supportRatio,"frequency ", max_freq, "center ", cent, "angle %.2f" % alfa, "axes (%.2f,%.2f)" % (a, bmax)
# removing founded ellipse pixels from further analysis
for p in supported:
pxy.remove(p)
# drawing founded ellipse
for i in range(EL_PATH_POINTS):
elx,ely = elData[i]
if i < EL_PATH_POINTS-1:
draw.line(elData[i] + elData[i+1], fill=colors[cIx])
io.save(sys.argv[2])
print "***************************************************************"
print "************************** DONE *******************************"
print "***************************************************************"

(Prototype algorithm is slow, tested only on 50x50 images). So, by running this algo on this image:
we will get such algorithm output:

Have fun in computer vision !

#### 1 comment:

1. Hello there,
I really like your blog. You have some really nice stuff here and I do hope you will carry on writing.
I am looking for passionate writers to join our community of bloggers - glipho.com
It might be a good idea to give your writing and your blog more exposure while having fun and meeting fellow writers.
Please check us out and drop me a line at hubert@glipho.com for any questions.
Best!
Hubert

Comment will be posted after comment moderation.