Who I'm awed by

Current favourites below. It is clearly a disservice to fit many of these people under a category, because clearly their influence and great talents are wide-ranging. In fact, the common theme in those that I'm awed by is as much breadth as depth. Salman Khan (of Khan Academy) it seems knows about everything well enough to teach it to you better than your specialist teacher. Bill Gates fluidly moves between high technology and philanthropy. Ramachandran is a brilliant and eccentric neuroscientist, but his ameteur attempts at archaeology are also very impressive. Tyler Cowen and Rajan look more like specialists in this group, but their breadth inside economics is wide ranging - unlike most economists who are limited to one or two areas. Some of it, I confess, is because of of my own inclinations in Economics. Daniel Kahneman is more of a specialist, yes, but in an area that has increased greatly our understanding of how breadth has a depth all its own. David Eagleman is an outstanding neuroscientist whose work of fiction "Sum" was equally outstanding. On a similar vein, Vikram Seth - delightful poet and mighty novelist; is also a pretty good painter; George Soros, who you might better know as a hedge fund manager, made seminal dents in the philosophy of markets. Charlie Rose is the most versatile interviewer I know - and betrays a good understanding of subjects from Literary theory to Macroeconomics to Hollywood to Brain sciences to Politics. In fact, he has been responsible for introducing me to several things I later got very interested in.

In time, I'm sure this list will change as my own interests evolve, but for now I'm pretty happy with this. There are so many others I admire a lot, but in condensing the list I know I've tended to favor those who are not just great, but whose greatness exerts itself across many disciplines, as well as into a human element I love.

Entrepreneurs:
Salman Khan
Bill Gates

Academics:
VS Ramachandran
Raghuram Rajan
Tyler Cowen
David Eagleman

Philosophy:
Iain McGilchrist
Jaggi Vasudev

Artists:
Vikram Seth
Jagjit Singh
Hugh Laurie

Eloquence:
Shashi Tharoor
Barack Obama

IOAM 2

I never quite know if I am replete or not, and it seems that most other people seem to know this. I do sense hunger. Not often - because I take my meals more or less 3 times a day so I end up not letting hunger a chance to set in - but on occasions when for whatever reasons I haven't been able to, I do feel, physically, pangs of hunger. It's not something to observe under IOAM because it is the most normal thing. I only note it because I've almost never quite known the counteracting feeling of experiencing, physically, that I am replete. I just stop out of conventional wisdom and reasonable rationing borne of external stimuli (e.g. gathered knowledge to the effect of "you're supposed to stop after eating this much.."). I wonder if it's just me, or most, or everyone.

IOAM 1

Coffee has never worked for me. Never made me more active, awake, anything. So far, it seems, probably never will. I like its taste though when mixed with milk and sugar.

Inchoate

The thoughts we entertain today enslave us tomorrow. If it is an identity, not a hypothesis, then the choice we have is not regarding whether we will be enslaved, but what we will be enslaved by.  

Calling a spade a spade

Rational discourse is standing on its last legs. After all, it is such a drag in front of impressive intellectual outbursts, especially when directed against somebody well-known. Calling oneself smart is cleverly understood to be in bad taste, but hey, there's an easy fix - calling others imbeciles is fine! It obviates the need to explicitly state how you're so much more evolved. Unfortunately, this quickest road to eternal smarts is gaining widespread traction. Who needs to do the hard work of understanding complicated policy implications, when calling Rahul Gandhi a pappu suffices, never mind your own secret failures at telling stare from stair.

 Let me state my disagreement with the extent of slamming Salman Khan, the bollywood actor, is facing for comparing his physically damaging routine while shooting for Sultan, a movie where he plays a wrestler, to rape. I'm no Salman Khan fan, but I'm a fan of trying to understand things in their context. Level-headedness and humility demand that one weigh their conclusions and metaphors before spewing them out to fawning followers. And that is also why I think the comparison Salman Khan drew was ill-advised. But the same preference for carefully weighing what you say is now surprisingly amiss, also, in those who castigate him in every possible way, for what was in all likelihood more a conversational gaffe than the all-out assault on women power it is now being portrayed to be.

 Intellectual heft, today, increasingly depends on identifying the prevailing fashions of the intelligentsia and magnifying that rhetoric in passably good prose. Thinking for oneself is a minority sport, and objectivity is for losers. By his admittedly imperfect usage, did Salman Khan really mean to "trivialize" rape, or did he mean to intensify his expression of the debilitating nature of his wrestling sequences? The answer to this question is not determined by morality or sentiment, but by plain old verbal comprehension. Yes, that boring thing that we only care about when preparing for CAT. Approaching it thus, Salman Khan used a simile - and similes don't have to be perfectly logical, or perfectly concordant. In a simile one often compares A to a more widely understood thing B, often to aggrandize the effect of A. In Hindi we refer to it as Atishyokti alankar - in English, I guess, one would call that hyperbole. As far as my comprehension serves me, in comparing his training/stunts to rape, he grossly over-aggrandized, and which I'm not defending. What I intend to bring forth, in fact, is that it was much more an aggrandizement (excessive, yes, and definitely unwarranted) of his stunts, than a downplaying or trivializing of rape. Efforts to understand and portray it as primarily the latter, in my opinion, distort the reality.

 Let me take the example of the talented singer Sona Mohapatra, of whom, unlike Salman Khan, I actually am a fan, and who took to social media recently to lambaste Salman Khan, calling him stupid, nasty, dangerous and an imbecile, for trivializing rape. Except that when you call a functioning individual an imbecile, you trivialize the hardships faced by the mentally challenged. "Ah, that's silly, you're twisting her intention." Well, maybe. So how about a song of hers she posted the very same day, titled qatl-e-aam, with its recurring "aap ki aankhon mein warna aaj qatl-e-aam hai". As heinous as rape obviously is, I'm sure qatl-e-aam, or mass murder, is just as bad. So is she trivializing murder when she compares it to some arbitrary dude's eyes, insensitive to all those whose loved ones have, in the long bloody history of time, been murdered?

 "Come on, now, you're making no sense - you're just pissed with Sona Mohapatra" you might be tempted to say at this point. And I would whole-heartedly agree - that the reasoning I provide can only be considered spite with Sona Mohapatra, not a cogent advocacy of victims of mass murders. I've already made up my mind on who I want to put down, and I'm stretching things way out of context to prove my point. Alas, that is precisely the point I've been trying to make.

 In a world full of politically-correct bandwagon-jumping, calling BS is a public good. Happy to do my part. Check her song out though, it is pretty awesome: https://www.youtube.com/watch?v=vosWN0CkDP4&sns=em

On cliches

The trouble with cliches is not their content, the trouble is the fact that they're cliched, so the message passes you by without registering itself. The trouble is that you fail to really take notice of it even once because it is floated around hundreds of times. In other words, the trouble with cliches is nothing, the trouble is with you. Notice.

Time Series in R .. to be continued

library(quantmod)

SBUX <- auto.assign="F)</p" from="2012-1-1" getsymbols="" to="2012-12-31">
plot(SBUX$SBUX.Adjusted, main="SBUX Closing Price (adjusted) for 2012")

// Typically there are many things we might want to do in time series analysis
// *generate forecasts of future values (predictive distributions)
// *gain an understanding of the process itself

//Constructing Indicators
//how can a ts of past info be "transformed" into indicators useful for forecasting

//Standard ts models fit a process to the series, requiring the estimation of a few parameters

//Standard approach
//1. Plot the time series. There are a lot of features that are still most readilty identified by a human
//2. Transform original ts to a stationary ts
//(a) determine whether to apply boxcox/log transformation or not
//(b) check for trend, seasonality, sharp changes, outliers, etc then do appropriate transformation
//3. Fit a model to the stationary ts
//4. Disgnostic tests to check model is reasonable, if its not go back to step 2
//5. Generate forecasts (predictive distributions) for the model
//6. Invert the transformation made in step 2, get forecasts of original ts.

# Stationarity

//Weakly stationary ts if :
//Var(X(t))//E[X(t)]=mu for all t ---------------------------------(2)
//autocov(r,r+h)=Cov(X(r), X(r+h))= f(h) = autocov(h) --(3)

//Also (3) implies that Var(X(t)) is independent of t, that is, a constant. So (1) can be restated.

//Examples
//(1) White noise process WN(mu, SigmaSq) is stationary
//(2) Random Walk process X(t) = Summ(1,t) [WN(0,SigmaSq)] is NOT stationary
//(3) MA(1) process X(t) = WN(t) + theta*WN(t-1) is stationary
//(4) ARCH(1) process is stationary, actually ARCH(1) process is also a white noise process


acf(SBUX$SBUX.Adjusted)

Box.test(SBUX$SBUX.Adjusted,lag=10,type="Ljung-Box")

pacf(SBUX$SBUX.Adjusted)

//-----------------------------

//AR 1

//Stationarity of AR1. Stationary if |fi1| < 1

//use arima() to fit

plot(SBUX$SBUX.Adjusted)
Box.test(SBUX$SBUX.Adjusted,lag=10,type="Ljung-Box", fitdf=1) //fitdf=1 since AR(1), prevents overfitting
arima(SBUX$SBUX.Adjusted, order=c(1,0,0))

plot(diff(log(SBUX$SBUX.Adjusted))[-1])
Box.test(diff(log(SBUX$SBUX.Adjusted))[-1],lag=10,type="Ljung-Box", fitdf=1)
arima(diff(log(SBUX$SBUX.Adjusted))[-1], order=c(1,0,0))

//AR(p)

DGS3 <- auto.assign="F)</p" getsymbols="" src="FRED">DGS3<-last p="" years="">head(DGS3diff)
DGS3diff<-diff 1="" p="">DGS3diff<-dgs3diff diff="" is.na="" p="">
pacf(DGS3diff) #AR1 not appropriate as in AR1 lags 2,3.. shouldn't matter after accounting for lag 1
#In general in PACF plot for AR(p), shouldn't see significant numbers after lag p --- Important

#acf formula is hard for AR(p) but there is function ARMAacf
#see how the ACF and PACF for an AR(2) model with fi1=0.5 and fi2=-0.2 looks:
plot(0:10, ARMAacf(ar=c(0.5,-0.2), lag.max=10), type="b", pch=16) #ACF
plot(1:10, ARMAacf(ar=c(0.5,-0.2), lag.max=10, pacf=TRUE), type="b", pch=16)  #PACF

# using auto.arima() to fit.
# 1. argument max.p sets the largest value of p to consider, max allowed is 5.
# 2. to fit a pure AR(p) model, either use arima() or in auto.arima() specify max.q=0
# 3. ic="aic" or "bic", or "aicc" which is corrected AIC - used for model selction by auto.arima()
# 4. input data to this fn should be vector, not xts. Apply as.numeric()
# 5. stepwise=FALSE to be specified if you want exhaustive search
# 6. seasonal=FALSE if don't want to consider models with seasonal components (default)
# 7. argument max.q sets the largest value of q to consider, max allowed is 5.
# 8. d can be specified explictly (e.g. d=2) else it will determine itself (not preferred).


# Consider daily trading volume of MSFT for 2012, find first difference of series.
# what order AR is chosen with aicc?

# Answer

library(forecast)
MSFT <- auto.assign="F)</p" from="2012-1-1" getsymbols="" to="2012-12-31">MSFTVolumeDiff <- diff="" olume="" p="">auto.arima(as.numeric(MSFTVolumeDiff), max.p=6, ic="aicc")

//It fit a model with p=1 and q=1. In addition, mean parameter was found to be unnecessary.

//Stationarity of AR(p)
//Recall that root of a fn f(x) is value of x that makes f(x) equal 0.
//It can be shown that AR(p) is stationary if EVERY root of f(x) = 1 - Summ:n(1,p)[FI(n)*x^n]
//has absolute value greater than 1.
//To get the roots of f(x) shown above we can use in R, polyroot(c(1, -FI(1), -FI(2),...))

# for example, lets say an AR(2) with FI(1) -1.2074 and FI(2) -0.2237 is fit. Is it stationary?
polyroot(c(1,-1.2074,0.2237))
# both are greater than 1, but the first one is pretty close.
# Ruppert says, "since the roots are estimated with error, reason to suspect nonstationarity"

# Unit Root Tests
# These are hypothesis tests for stationarity
# the simplest one being Dickey-Fuller test - but useful only for AR(1)
# Augmented Dickey Fuller (ADF) is more commonly used: adf.test()

library(tseries)
adf.test(as.numeric(DGS3diff), alternative="s") #"s" for stationarity, "e" for explosive
#p-value of 0.01 suggests we can reject the null in favor of alternative, that is "stationary"

# Lets see what auto.arima fits on DGS3diff
auto.arima(as.numeric(DGS3diff), max.p=6, ic="aicc")
# It fit ARIMA(2,0,2) with some parameters. Can't use polyroot on it since it has MA component.

// MA, ARMA, ARIMA

# MA(q) process - Stationary, no correlation beyond lag q in ACF plot.

# Lets simulate MA(2) process data of length 250 with theta1=0.3 and theta2=-0.2 and sigmaSq=1
simdat <- arima.sim="" model="list(ma=c(0.3,-0.2))," n="250," sd="sqrt(1))</p">
# See if the simulated data is stationary by using ADF test
adf.test(as.numeric(simdat),alternative="s") #small p value suggests stationary data

# See if any correlation beyond lag 2 in ACF plot
acf(simdat) # no correlation beyond lag 2, as expected, since data was simulated from MA(2) process

# Let's fit a MA(2) model to this data (as if we didn't know it was simulated data from specified MA(2))
holdfit <- arima="" order="c(0,0,2))</p" simdat="">holdfit # returns an MA(2) with theta1 0.24, theta2 -0.16 and sigmaSq 0.88

# But the theta1 estimate has a std error - lets see the confidence interval for it
holdfit$coef[1] - 1.96*sqrt(holdfit$var.coef[1,1])
holdfit$coef[1] + 1.96*sqrt(holdfit$var.coef[1,1]) # 95% C.I for theta1 is (0.11,0.35)

#Ljung-Box Test (H1: "at least some correl at some lag") applied to residuals of the fit
Box.test(holdfit$residuals, fitdf=2, lag=20, type="Ljung-Box") #p not small enough to reject null-> no correl
acf(holdfit$residuals) # acf confirms box test result
qqnorm(holdfit$residuals) # normal probability plot.. need to make this work.
qqline(holdfit$residuals)

# Comparing acf from simulated data to acf from true process
acf(simdat, cex.axis=1.3,cex.lab=1.3,main="") # Sample ACF of simulated data
trueacf=ARMAacf(ma=c(0.3,-0.2), lag.max=30) # true ACF
points(0:30,trueacf, pch=16, col=2)
fitacf=ARMAacf(ma=holdfit$coef[1:2],lag.max=30) #ACF of the fitted model
points(0:30,fitacf, col=4)
legend(23,1,legend=c("True ACF","Fitted model ACF"), pch=c(16,1),col=c(2,4), xjust=1, cex=1.3)

# fitting ARIMA

DTWEXM <- auto.assign="F)" dollar="" getsymbols="" index="" p="" src="FRED" trade-weighted="">DTWEXM <- dtwexm="" p="">plot(DTWEXM) # clearly not stationary
plot(diff(DTWEXM,1,1)) #diff(vec,n,d) is lag n difference (X(t)-X(t-n)) applied d times #now looks stationary
plot(diff(DTWEXM,1,2)) #2nd (lag-1) differences .. this too looks stationary.

# We want smallest d that's good enough
adf.test(as.numeric(diff(DTWEXM,1,1)[!is.na(diff(DTWEXM,1,1))]), alternative="s")
adf.test(as.numeric(diff(DTWEXM,1,2)[!is.na(diff(DTWEXM,1,2))]), alternative="s")
# Looks like 1st differencing was sufficient, so d=1

# REVISION: using auto.arima() to fit.
# 1. argument max.p sets the largest value of p to consider, max allowed is 5.
# 2. to fit a pure AR(p) model, either use arima() or in auto.arima() specify max.q=0
# 3. ic="aic" or "bic", or "aicc" which is corrected AIC - used for model selction by auto.arima()
# 4. input data to this fn should be vector, not xts. Apply as.numeric()
# 5. stepwise=FALSE to be specified if you want exhaustive search
# 6. seasonal=FALSE if don't want to consider models with seasonal components (default)
# 7. argument max.q sets the largest value of q to consider, max allowed is 5.
# 8. d can be specified explictly (e.g. d=2) else it will determine itself (not preferred).
# 9. approx=F tso that true maximum likelihood estimates are found

fitted_model<-auto .arima="" a="" arima="" as.numeric="" d="1," fit="" ic="aic" it.="" it="" p="" see="" that="" to="" we="">
#---------------------------------------------------------------------------------------

# Forecasting

# REVISION: Standard approach
//1. Plot the time series. There are a lot of features that are still most readilty identified by a human
//2. Transform original ts to a stationary ts
//(a) determine whether to apply boxcox/log transformation or not
//(b) check for trend, seasonality, sharp changes, outliers, etc then do appropriate transformation
//3. Fit a model to the stationary ts
//4. Disgnostic tests to check model is reasonable, if its not go back to step 2
//5. Generate forecasts (predictive distributions) for the model
//6. Invert the transformation made in step 2, get forecasts of original ts.

k=3 # number of future predictions
# Rem that holdfit was the model fitted to MA(2) simulated values
next_k_preds <-predict holdfit="" n.ahead="k)</p">
# Lets do another example

AAA <- auto.assign="F," getsymbols="" p="" src="FRED">AAA <- aaa="" p="">
plot(AAA) # Clearly not stationary, there is a trend
plot(diff(AAA)) # looks like trend has been removed, could be stationary
adf.test(as.numeric(diff(AAA)[-1]),alt="s") # small p value suggests stationary
holdfit <- a="" approx="F)" arima="" as.numeric="" auto.arima="" d="1," fits="" it="" p="" process="">
#diagnostic tests
Box.test(holdfit$residuals, fitdf=1, lag=10, type="Ljung-Box") #high p-value, fail to reject null: no correl
acf(holdfit$residuals)# suggests no serial correlation among residuals
pacf(holdfit$residuals)# suggests no serial conditional correlation among residuals
plot(holdfit$residuals) # looks good

#generate forecasts
holdpreds <- 6="" also="" care="" holdfit="" inverting="" n.ahead="6)" note="" of="" p="" predict="" r="" step="" takes="" that="" the="" transformation="">holdpreds

#--------------------------------------------------------------------------------------------

# Seasonality

HOUSTNSA <- auto.assign="F)</p" getsymbols="" src="FRED">HOUSTNSA <- houstnsa="" p="">plot(HOUSTNSA, ylab="Housing Starts (thousands of units)")

# We notice that there are trends not stationary
# We also notice seasonality
acf(HOUSTNSA) # confirms seasonality, peaks every 12 lags

plot(diff(HOUSTNSA)[-1]) #this removes the trend, but seasonality remains acf(diff(HOUSTNSA)[-1])
plot(diff(HOUSTNSA, 12,1)) #this removes seasonality but trends remain. (seasonally adjusted)
plot(diff(diff(HOUSTNSA, 12,1))) # this removes both and should be good. ------------(A)
# we can fit ARMA to it, but in practice this turns out to not yield models with sufficient flexibility

# Multiplicative ARIMA models for accomodating seasonality (preferred)
# ARIMA((p,d,q)*(P,D,Q)s) where s is the period of seasonality (12 for a year, for example)
# Think of it as fitting ARIMA(p,d,q) on non-seasonal component and ARIMA(P,D,Q) on seasonal component of ts.
# This multiplicative ARIMA model is fitted upon the ts obtained in (A) above
#(In R you needn't do the differencings first)

# Fitting a multiplicative ARIMA model

# METHOD 1: using arima() fitting a ARIMA((1,1,1)*(1,1,0)12) to the housing data again.
holdfit1 <- arima="" order="c(1,1,1)," period="12))</p" seasonal="list(order=c(1,1,0),">
# METHOD 1: using auto.arima()

# s is typically already known, or estimated exploring the ACF, not something auto.arima() finds.
# ts object fed to auto.arima should already be in a form that specifies the period of the series, i.e. s

HOUSts<-ts an="" creates="" frequency="12)" function="" not="" object="" p="" r="" series="" the="" time="" ts="" xts=""># d, D should always be specified yourself based on inspection of series, although auto.arima can search for it
holdfit2<-auto .arima="" d="1,D=1,approx=F)</p" ts="">
# We got a better (smaller) aic than that in the totally self-specified arima() model holdfit1.

# Use it for forecasting now.
holdpreds <- holdfit2="" n.ahead="24)</p" predict="">plot(HOUSts, xlab="Year", xlim=c(0,22))
lines(holdpreds$pred,col=2)
lines(holdpreds$pred + (1.96*holdpreds$se), col=2, lty=2)
lines(holdpreds$pred - (1.96*holdpreds$se), col=2, lty=2)

#---------------------------------------------------------------------------------------------

# Heteroskedasticity

#---------------------------------------------------------------------------------------------

# ARFIMA models

#---------------------------------------------------------------------------------------------

# Revisiting Regression

# One key assumption in regression was that errors iid normal, i.e, uncorrelated with const variance
holdlm=lm(HOUSts~sqrt(HOUSts)) # just some random regression

# In a time series regression test whether they are really uncorrelated makes sense to see ACF of residuals
acf(holdlm$residuals, lag.max=10)

# Another way is the durbin-watson hypothesis test, made specifically to test serial correl in a linear model
# H0: there is no correlation. If small p-value, reject the null: there is correlation.
library(car)
durbinWatsonTest(holdlm,10)

# DW suffers from the same multiple testing problem as acf. Can use Ljung-Box on holdlm$residuals, if needed.
# In our example we can see that errors ARE correlated.

# How to deal with this problem?
# 1. might want to log transform the predictor or response variable or both in the linear model
# 2. If significant serial correlation, fit an ARIMA(p,d,q) process in place of model error, instead of WN.

# Example. Response var is 3-year tsy rate, predictor var is 1-yr tsy-rate.
OneYearTreas <- auto.assign="F)</p" getsymbols="" src="FRED">ThreeYearTreas<-getsymbols auto.assign="F)</p" src="FRED">
holdfit1 <- hreeyeartreas="" lm="" oneyeartreas="" p="">holdfit2 <- arima="" hreeyeartreas="" order="c(1,1,1)," xreg="OneYearTreas)</p">holdfit3 <- as.numeric="" auto.arima="" d="1)</p" hreeyeartreas="" xreg="as.numeric(OneYearTreas),">acf(holdfit3$residuals)

# This solves the problem of serial correlation. What about non-constant variance?
# One way is model the errors as GARCH. Unfortunately not straightforward like with ARIMA above. Follow steps:
#1. Fit a regular lm model.
#2. Fit a time series model to the residuals  from model fit above.
#3. Refit the model using weighted least squares, weights being reciprocals of conditional variances from step2
#   Conditional variance from garchFit() are found as: holdgarchfit@sigma.t
#   weights = 1/(holdfit@sigma.t^2)
# Do it yourself.

#----------------------------------------------------------------------------------------------