Hey Calvin, ACFs produced by your code are biased towards zero...
The reason for that is that the first k elements in cross_sum (variable of the sub-function) are always zero. Also, dimensions of cross_sum after the loop in lines 104-106 are always Nx1. In large sample the bias is small but in small samples it might be sensible.
Given that matlab is very bad at handling loops it is better to avoid them altogether if possible. I adjusted your code by removing the sub-function completely, "global" attributes for N and ybar (lines 46 and 48) and substituting loop in lines 52-54 by
for i = 1:p
yvar = (y-ybar)'*(y-ybar) ;
ta(i) = (cross_sum / yvar)*(N/(N-i)) ;