We introduce a novel non-parametric methodology to test for the dynamical time evolution of the lag-lead structure between two arbitrary time series. The method consists of constructing a distance matrix based on the matching of all sample data pairs between the two time series. Then, the lag-lead structure is searched for as the optimal path in the distance matrix landscape that minimizes the total mismatch between the two time series, and that obeys a one-to-one causal matching condition. To make the solution robust to the presence of a large amount of noise that may lead to spurious structures in the distance matrix landscape, we generalize this optimal search by introducing a fuzzy search by sampling over all possible paths, each path being weighted according to a multinomial logit or equivalently Boltzmann factor proportional to the exponential of the global mismatch of this path. We present the efficient transfer matrix method that solves the problem and test it on simple synthetic examples to demonstrate its properties and usefulness compared with the standard running-time cross-correlation method. We then apply our 'optimal thermal causal path' method to the question of the lag-dependence between the US stock market and the treasury bond yields and confirm our earlier results on an arrow of the stock markets preceding the Federal Reserve Funds' adjustments, as well as the yield rates at short maturities in the period 2000-2003. Our application of this technique to inflation, inflation change, GDP growth rate and unemployment rate unearths non-trivial lag relationships: the GDP changes lead inflation especially since the 1980s, inflation changes leads GDP only in the 1980 decade, and inflation leads unemployment rates since the 1970s. In addition, our approach seems to detect multiple competing lag structures in which one can have inflation leading GDP with a certain lag time and GDP feeding back/leading inflation with another lag time.