# Time Series extension to APPL (http://applsoftware.com)
# Keith H. Webb
# khwebb@math.wm.edu

# Express an ARMA(p, q) model as Y[t] = c + phi[1]y[t-1] + ...
ARMAModel := proc(p :: nonnegint, q :: nonnegint, lag :: nonnegint)
local y, i:
if (nargs <> 3) then
  print(`ERROR(ARMAModel): This procedure requires 3 arguments:`):
  print(`The parameters p, q in ARMA(p, q)`):
  print(`The lag, i.e. k such that the equation is y[t-k]`):
  RETURN():
fi:

y[t-lag] := c:
if (p > 0) then
  for i from 1 to p do
    y[t-lag] := y[t-lag] + phi[i]*y[t-i-lag]:
  od:
fi:

y[t-lag] := y[t-lag] + epsilon[t-lag]:
if (q > 0) then
  for i from 1 to q do
    y[t-lag] := y[t-lag] + theta[i]*epsilon[t-i-lag]:
  od:
fi:

RETURN(y[t-lag]):
end:

# Take an APPL distribution in string form and an epsilon[t-j] term and find the new
#  distribution of epsilon[t-j] at time t-j.
with(StringTools):
NewDist := proc(term, Stringdist)
local temp, temp2, Y:
global t:

temp := Substitute(term, "^2", ""):
temp := Substitute(temp, "epsilon[", ""):
temp := Substitute(temp, "t", ""):
temp := Substitute(temp, "-", ""):
temp := Substitute(temp, "]", ""):        # isolate j
if (temp = "") then
  temp := 0:
else
  temp := parse(temp):
fi:
temp2 := convert(t - temp, 'string'):
temp2 := Insert(temp2, 0, "("):
temp2 := Insert(temp2, length(temp2), ")"):
Y := SubstituteAll(Stringdist, "tau", temp2):     # replace all taus in the APPL dist. with t-j
Y := SubstituteAll(Y, "~", ""):                   # get rid of ~, Maple's way of marking variables with assumptions
assume(t > temp):
Y := parse(Y):
RETURN(Y):
end:

# Shorthand procedures for creating the data structure
MA := proc(q :: {posint, list(realcons)})
RETURN(ARMA(0, args)):
end:

AR := proc(p :: {posint, list(realcons)})
if (nargs = 2) then
  RETURN(ARMA(p, 0, args[2])):
else
  RETURN(ARMA(p, 0)):
fi:
end:

ARMA := proc(p :: {nonnegint, list(realcons)}, q :: {nonnegint, list(realcons)})
local X:

if (p = 0 and q = 0) then
  print(`ERROR(ARMA): At least one of p or q must be positive.`):
  RETURN():
fi:

X := [p, q]:
if (nargs = 3) then
  X := [op(X), args[3]]:
fi:
X := TSParse(X):

if (type(p, list(realcons))) then
  X[1] := p:
fi:
if (type(q, list(realcons))) then
  X[2] := q:
fi:

RETURN(X[1..3]):
end:

# Take the data structure used by users and convert it into the form used internally, making
#  the first two arguments integers rather than lists and adding "time-dependent" or "time-independent"
TSParse := proc(X)
local r, temp:

r := [0, 0]:
if (type(X[1], nonnegint) and type(X[2], nonnegint)) then
  r := [X[1], X[2]]:
elif (type(X[1], list(realcons)) and type(X[2], list(realcons))) then
  r := [nops(X[1]), nops(X[2])]:
elif (((X[1] = 0) xor (X[2] = 0)) and (type(X[1], list(realcons)) xor type(X[2], list(realcons)))) then
  if (X[1] = 0) then
    r := [0, nops(X[2])]:
  else
    r := [nops(X[1]), 0]:
  fi:
else
  print(`ERROR: Time series must be given in a list format, etc.`):
  RETURN():
fi:

if (nops(X) > 2 and X[3] <> 0) then
  r := [op(r), X[3]]:
  temp := convert(X[3], 'string'):
  if (Search("tau", temp) = 0) then
    r := [op(r), "time-independent"]:
  else
    r := [op(r), "time-dependent"]:
  fi:
else
  print(`Error terms are assumed to be N(0, 1)`):
  r := [op(r), NormalRV(0, 1), "time-independent"]:
fi:
RETURN(r):
end:

# Find the mean of a time series. This differs from c when the error term distribution has a
#  non-zero mean. Essential for finding variance and co-variance.
TSMean := proc(TS :: list, lag :: nonnegint)
local xt, X, Y, tsmean, ystring, ylist, Emean, termmean, i, j, temp, termlist, tau, time:

if (nargs <> 2 and nargs <> 3) then
  print(`ERROR(TSMean): This procedure requires 2 arguments:`):
  print(`A time series model written as a sum.`):
  print(`The distribution of the error terms (epsilon terms)`):
  RETURN():
fi:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (xt[1] > 0) then
  if (nargs > 2 and type(args[3], posint)) then
    time := args[3]:
  else
    time := readstat("Find mean for y[t], t="):
  fi:
  ystring := convert(ARReduce(xt[1], xt[2], time), 'string'):
else
  ystring := convert(ARMAModel(xt[1], xt[2], lag), 'string'):
fi:

X := convert(xt[3], 'string'):
tsmean := 0:
ylist := Split(ystring, "+"):

# The idea of this procedure is to the time series as a string into a list at each +, then
#  split all terms by * and find expected value on a term-by-term basis
for i from 1 to nops(ylist) do
  termmean := 1:
  termlist := Split(ylist[i], "*"):
  for j from 1 to nops(termlist) do
    if (IsSubSequence("epsilon", termlist[j])) then
      if (xt[4] = "time-dependent") then
	Y := NewDist(termlist[j], X):
      else
	Y := xt[3]:
      fi:
      termmean := termmean * Mean(Y):
    else
      temp := parse(termlist[j]):
      termmean := termmean * temp:
    fi:
  od:
  tsmean := tsmean + termmean:
od:

RETURN(TSSubstitute(TS, tsmean)):
end:

# Find the variance of a time series
TSVariance := proc(TS :: list, lag :: {nonnegint, list(nonnegint)})
local xt, ystring, ylist, i, yvar, termvar, temp, X, Y, U, j, g, yexpand, termlist, xmean, ymean, xtemp, ytemp, lags, time:

if (nargs <> 2 and nargs <> 3 and nargs <> 4) then
  print(`ERROR(TSVariance): This procedure requires 2 arguments:`):
  print(`An ARMA model in list form`):
  print(`The lag of the autocovariance`):
  RETURN():
fi:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (member("YW", [args]) and xt[1] > 0) then
  RETURN(YuleWalker(TS, lag, "var")):
fi:

if (nops(lag) < 2) then # lag is a list in the case where the process is not covariance stationary
  lags := [0, lag]:     #  and variance between y[t-j] and y[t-j] needs to be found
else
  lags := lag:
fi:

if (xt[1] > 0) then
  if (nargs > 2 and type(args[3], posint)) then
    time := args[3]:
  else
    time := readstat("Find (co)variance for y[t], t="):
  fi:
  xtemp := ARReduce(xt[1], xt[2], time - lags[1]):
  ytemp := ARReduce(xt[1], xt[2], time - lags[2]):
else
  xtemp := ARMAModel(xt[1], xt[2], lags[1]):
  ytemp := ARMAModel(xt[1], xt[2], lags[2]):
  time := 0:
fi:

xmean := TSMean(xt, lags[1], time):
ymean := TSMean(xt, lags[2], time):
yexpand := expand((xtemp - xmean) * (ytemp - ymean)):
ystring := convert(yexpand, 'string'):
ystring := Insert(ystring, 0, "0+"):
ystring := SubstituteAll(ystring, "-", "+(-1)*"):
ystring := SubstituteAll(ystring, "t+(-1)*", "t-"):
ystring := SubstituteAll(ystring, "++", "+"):
ystring := SubstituteAll(ystring, "/", "*1/"):
ystring := RegSubs("/\\(\\+\\(-1+\\)\\*" = "/\(-", ystring):
ystring := RegSubs("(\\([-0-9]+)\\+([0-9\\*]*t\\))" = "\\1-\(-\\2\)", ystring):    # Yuck!!!
ystring := RegSubs("(\\([-0-9*]*t)\\+([0-9]+\\))" = "\\1-\(-\\2\)", ystring):
ylist := Split(ystring, "+"):

X := convert(xt[3], 'string'):
yvar := 0:

# As in TSMean, the expaned string is split so that expected value can be computed term-by-term
for i from 1 to nops(ylist) do
  termvar := 1:
  termlist := Split(ylist[i], "*"):
  for j from 1 to nops(termlist) do
    if (not IsPrefix("epsilon", termlist[j])) then # Parse constant terms
      temp := parse(termlist[j]):
      termvar := termvar * temp:
    fi:
    if (IsPrefix("epsilon", termlist[j]) and IsSuffix("^2", termlist[j])) then # find mean of epsilon[t]^2
      if (xt[4] = "time-dependent") then
	Y := NewDist(termlist[j], X):
      else
	Y:= xt[3]:
      fi:
      g := [[x -> x^2, x -> x^2], [-infinity, 0, infinity]]:
      U := Transform(Y, g):
      termvar := termvar * Mean(U):
    fi:
    if (IsPrefix("epsilon", termlist[j]) and not IsSuffix("^2", termlist[j])) then # find mean of epsilon[t]
      if (xt[4] = "time-dependent") then
	Y := NewDist(termlist[j], X):
      else
	Y:= xt[3]:
      fi:
      temp := Mean(Y):
      termvar := termvar * temp:
    fi:
  od:
  yvar := yvar + termvar:
od:
RETURN(TSSubstitute(TS, yvar)):
end:

# Find autocorrelation, "YW" flag optional
TSCorrelation := proc(TS :: list, lag :: nonnegint)
local xt, devx, devy, cov, temp, corr, time:

if (nargs <> 2 and nargs <> 3 and nargs <> 4) then
  print(`ERROR(TSCorrelation): This procedure requires 3 arguments:`):
  print(`Two time series models written as sums.`):
  print(`The distribution of the error terms (epsilon terms)`):
  RETURN():
fi:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (lag = 0) then
  RETURN(1):
fi:

if (member("YW", [args]) and xt[1] > 0) then
  RETURN(YuleWalker(TS, lag)):
fi:
if (xt[1] > 0) then
  if (nargs > 2 and type(args[3], posint)) then
    time := args[3]:
  else
    time := readstat("Find correlation for y[t], t="):
  fi:
else
  time := 0:
fi:

cov := TSVariance(xt, lag, time):
devx := sqrt(TSVariance(xt, [0, 0], time)):
devy := sqrt(TSVariance(xt, [lag, lag], time)):
temp := cov / (devx * devy):
corr := simplify(temp):
RETURN(TSSubstitute(TS, corr)):
end:

# Find partial autocorrelation
with(LinearAlgebra):
TSPartialCorrelation := proc(TS :: list, lag :: nonnegint)
local xt, time, rho, M, N, i, j, pacf, YW:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

YW := "":
if (member("YW", [args])) then
  YW := "YW":
fi:
if (xt[1] > 0) then
  if (nargs > 2 and type(args[3], posint)) then
    time := args[3]:
  elif (YW <> "YW") then
    time := readstat("Find partial correlation for y[t], t="):
  fi:
else
  time := 0:
fi:

if (lag < 2) then
  RETURN(TSCorrelation(TS, lag, time, YW)):
fi:

rho := [seq(TSCorrelation(TS, k, time, YW), k=0..lag)]:

M := Matrix(lag):
N := Matrix(lag):
for i from 1 to lag do:
  for j from 1 to lag do:
    M[i,j] := rho[modp(max(i,j) - min(i,j), lag) + 1]:
    N[i,j] := rho[modp(max(i,j) - min(i,j), lag) + 1]:
  od:
od:
for i from 1 to lag do:
  N[i,lag] := rho[i + 1]:
od:

pacf := Determinant(N) / Determinant(M):
RETURN(TSSubstitute(TS, pacf)):
end:

# If a result is given with aribtrary values and the inputted data structure has real parameters,
#  then this procedure performs the substitution
TSSubstitute := proc(TS, corr)
local i, temp, temp2, stringcorr:

stringcorr := convert(corr, 'string'):
if (type(TS[1], list(realcons))) then
  for i from 1 to nops(TS[1]) do
    temp := convert(TS[1][i], 'string'):
    temp := cat("(", temp, ")"):
    temp2 := cat("phi[", convert(i, 'string'), "]"):
    stringcorr := SubstituteAll(stringcorr, temp2, temp):
  od:
fi:
if (type(TS[2], list(realcons))) then
  for i from 1 to nops(TS[2]) do
    temp := convert(TS[2][i], 'string'):
    temp := cat("(", temp, ")"):
    temp2 := cat("theta[", convert(i, 'string'), "]"):
    stringcorr := SubstituteAll(stringcorr, temp2, temp):
  od:
fi:
RETURN(parse(stringcorr)):
end:

# Plots autocorrelations, even time-dependent ones.
with(plots):
TSPlot := proc(TS :: list)
local time, time2, graph, i, lags, points, theta, pstring, temp, temp2, xt, partial, lab, YW:
global c:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

c := 0: #AR finite-time autocorrelations often involve c, can't plot expression with variables
YW := "":
lags := 2*max(xt[1], xt[2]):
time := lags + 1:
time2 := time:
partial := 0:
if (nargs = 1) then
  if (xt[4] = "time-dependent") then
    print(`t is assumed to be`, time):
  fi:
elif (type(args[2], 'posint')) then
  if (xt[4] = "time-dependent") then
    time2 := args[2]:
  else
    time := args[2]:
    time2 := time:
  fi:
elif (args[2] <> "partial" and args[2] <> "YW") then
  print(args[2], `must be of type 'posint'`):
  RETURN():
fi:
if (member("partial", [args])) then
  partial := 1:
fi:
if (member("YW", [args])) then
  YW := "YW":
fi:

temp := 0:
if (time <> time2) then
  temp := 1
fi:

if (partial = 0) then
  lags := max(xt[1], xt[2]) + 3*xt[1]:
  points := [seq([k,TSCorrelation(TS, k, time, YW)],k=temp..lags)]:
else
  if (xt[2] <> 0) then
    lags := max(xt[1], xt[2]) + 3*xt[2]:
  else
    lags := xt[1]:
  fi:
  points := [seq([k,TSPartialCorrelation(TS, k, time, YW)],k=temp..lags)]:
fi:

lab := ["Lag", "Correlation"]:
if (partial = 1) then
  lab[2] := "Partial Correlation":
fi:
if (temp = 1) then
  lab[1] := "Time":
  for i from 1 to lags do
    points[i][1] := t:
  od:
fi:

pstring := convert(points, 'string'):
#for i from 1 to xt[1] do
#  if (TS[1][i] < 0) then
#    temp2 := convert(-TS[1][i], 'string'):    Changed xt to TS in the above calls to TSCorrelation and
#    temp := cat("(-", temp2, ")"):             TSPartialCorrelation so actual values are found there.
#  else                                         The string was too long for Maple to manipulate and was
#    temp := convert(TS[1][i], 'string'):       causeing errors.
#    temp := cat("(", temp, ")"):
#  fi:
#  pstring := SubstituteAll(pstring, cat("phi[", i, "]"), temp):
#od:
#for i from 1 to xt[2] do
#  if (TS[2][i] < 0) then
#    temp2 := convert(-TS[2][i], 'string'):
#    temp := cat("(-", temp2, ")"):
#  else
#    temp := convert(TS[2][i], 'string'):
#    temp := cat("(", temp, ")"):
#  fi:
#  pstring := SubstituteAll(pstring, cat("theta[", i, "]"), temp):
#od:

graph := [seq([0, 0], i=1..time2)]:
for i from time to time2 do
  if (time <> time2) then
    temp := SubstituteAll(pstring, "t", convert(i, 'string')):
    temp2 := parse(temp):
  else
    if (xt[4] = "time-dependent") then
      temp := SubstituteAll(pstring, "t", convert(i, 'string')):
      points := parse(temp):
    fi:
    RETURN(pointplot(points, labels=lab, axes=FRAME, xtickmarks=[seq(k, k=0..lags)])):
  fi:
  graph[i] := pointplot(temp2):
od:
graph := op(time..time2, graph):

RETURN(display(graph, labels=lab, axes=FRAME, xtickmarks=[seq(k, k=time..time2)])):
end:

#Performs unit root analysis
UnitRoots := proc(TS :: list)
local eqn, neweqn, newmodel, R, RA, RB, A, B, i, badroot, r, xt, temp:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (xt[4] = "time-dependent") then
  print(`ERROR(UnitRoots): Error terms must be time independent`):
  RETURN():
fi:
if (not type(TS[1], list(realcons)) and (TS[1] <> 0)) then
  print(`ERROR(UnitRoots): Parameters must be in a list of real values`):
  RETURN():
fi:
if (not type(TS[2], list(realcons)) and (TS[2] <> 0)) then
  print(`ERROR(UnitRoots): Parameters must be in a list of real values`):
  RETURN():
fi:

eqn := 1:
# Solve AR characteristic equation if model has AR parts
if (xt[1] <> 0) then
  for i from 1 to xt[1] do
    eqn := eqn - TS[1][i]*z^i:
  od:
  if (xt[1] > 4) then
    R := [evalf(solve(eqn = 0))]:
  else
    R := [solve(eqn = 0)]:
  fi:
  R := sort(R, (w, z) -> abs(evalf(w)) < abs(evalf(z))):
  RA := [seq(abs(evalf(R[i])), i=1..xt[1])]:
  RB := [seq([Re(evalf(R[i])), Im(evalf(R[i]))], i=1..xt[1])]:
  print(`The roots are:`):
  print(R):
  if (RA[1] > 1) then
    print(`All roots lie outside the unit circle, so the ARMA process is stationary`):
  else
    print(`The ARMA process is non-stationary`):
    badroot := xt[1]:
    for i from 1 to xt[1] do
      if (RA[i] > 1) then
	badroot := i - 1:
	temp := cat(`Number of roots inside the unit circle: `, convert(badroot, 'string')):
	print(temp):
	break:
      fi:
    od:
  fi:
  if (badroot = xt[1]) then
    print(`All roots lie inside the unit circle`):
  fi:
  A := pointplot(RB, symbol=circle):
  B := polarplot([1, t, t=0..2*Pi], color=blue):
  print(display({A, B}, labels=["R", "i"], scaling=CONSTRAINED)):
  RETURN():
fi:

for i from 1 to xt[2] do
  eqn := eqn + TS[2][i]*z^i:
od:

if (xt[2] > 4) then
  R := [evalf(solve(eqn = 0))]:
else
  R := [solve(eqn = 0)]:
fi:

R := sort(R, (w, z) -> abs(evalf(w)) < abs(evalf(z))):
RA := [seq(abs(evalf(R[i])), i=1..xt[2])]:
RB := [seq([Re(evalf(R[i])), Im(evalf(R[i]))], i=1..xt[2])]:

print(`The roots are:`):
print(R):

r := TS[2]:
if (RA[1] > 1) then
  badroot := 0:
  neweqn := eqn:
  print(`All roots lie outside the unit circle, so the MA process is invertible.`):
else
  badroot := xt[2]:
  for i from 1 to xt[2] do
    if (RA[i] > 1) then
      badroot := i - 1:
      temp := cat(`Number of roots inside the unit circle: `, convert(badroot, 'string')):
      print(temp):
      break:
    fi:
  od:
  neweqn := 1:
  for i from 1 to badroot do
    neweqn := neweqn * (1 - R[i]*z):
  od:
  if (badroot < xt[2]) then
    for i from (badroot + 1) to xt[2] do
      neweqn := neweqn * (1 - z/R[i]):
    od:
  else
    print(`All roots are inside the unit circle.`):
  fi:
  neweqn := expand(neweqn):
  newmodel := c + epsilon[t]:
  for i from 1 to xt[2] do
    r[i] := Re(coeff(neweqn, z, i)):
    if (abs(evalf(r[i])) < .00001) then
      r[i] := 0:
    fi:
    newmodel := newmodel + r[i]*epsilon[t-i]:
  od:
  print(`Assuming error terms are i.i.d., then an invertible representation is:`):
  temp := Y[t] = newmodel:
  print(temp):
fi:

A := pointplot(RB, symbol=circle):
B := polarplot([1, t, t=0..2*Pi], color=blue):
print(display({A, B}, labels=["R", "i"], scaling=CONSTRAINED)):

RETURN(r):
end:

# Performs the recursive substitution to transform, for example, an AR(8) model to an expression
#  involving only epsilons and phis.
ARReduce := proc(p, q, time)
local i, temp, temp2, TS, newp, newq:
global t:

newp := p:
newq := q:
t := time:
temp := ARMAModel(p, q, 0):
TS := convert(temp, 'string'):

for i from 1 to time do
  if (newp+i > time) then
    newp := newp - 1:
  fi:
  if (newq+i > time) then
    newq := newq - 1:
  fi:
  temp := ARMAModel(newp, newq, i):
  temp := cat("(", convert(temp, 'string'), ")"):
  temp2 := convert(y[time-i], 'string'):
  TS := SubstituteAll(TS, temp2, temp):
od:

unassign('t'):
RETURN(expand(parse(TS))):
end:

# Realizes an ARMA model with normal error terms via the Box-Muller transform
TSRealization := proc(TS :: list)
local xt, a, b, z, z1, z2, epsilon, eterms, i, j, n, theta, y, temp, iterations:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (nargs = 1) then
  iterations := 50:
elif (type(args[2], posint)) then
  iterations := args[2]:
else
  print(`ERROR(TSRealization): the second argument must be of type 'posint'`):
fi:

if (xt[4] = "time-dependent") then
  print(`Realization for time-dependent distributions are not supported yet`):
  RETURN():
fi:

theta := [1, op(TS[2])]:
n := iterations + max(2*xt[1], xt[2]):
z := array(1..n):
for i from 1 to ceil(n / 2) do
  a := rand()/10^12:
  b := rand()/10^12:
  z1 := sqrt(-2*ln(a)) * cos(2*Pi*b):      # Box-Muller
  z2 := sqrt(-2*ln(a)) * sin(2*Pi*b):
  z[i] := evalf(z1):
  z[n+1-i] := evalf(z2):
od:
z := convert(z, list):

eterms := Mean(xt[3]) + sqrt(Variance(xt[3])) * z:
y := array(1..n):
if (xt[1] > 0) then
  for i from 1 to max(xt[1], xt[2]) do
    y[i] := eterms[i]:
  od:
fi:
for i from max(xt[1]+1, xt[2]+1) to n do
  y[i] := 0:
  if (xt[1] > 0) then
    for j from 1 to xt[1] do
      y[i] := y[i] + TS[1][j]*y[i-j]:
    od:
  fi:
  for j from 0 to xt[2] do
    y[i] := y[i] + theta[j+1]*eterms[i-j]:
  od:
od:
y := convert(y, list):
print(listplot(y[n-iterations+1..n], axes=FRAME)):

RETURN():
end:

# Creates a spectral density function
SDF := proc(TS :: list)
local i, xt, f, time, YW, temp:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

YW := "":
if (member("YW", [args]) and xt[1] > 0) then
  YW := "YW":
fi:

if ((xt[1] > 0) and YW = "") then
  if (nargs > 1 and type(args[2], posint)) then
    time := args[2]:
  else
    time := readstat("Find SDF for y[t], t="):
  fi:
else
  time := 0:
fi:

f := TSVariance(TS, 0, time, YW):
for i from 1 to max(1000*xt[1], xt[2]) do
  temp := TSVariance(TS, i, time, YW):
  f := f + 2*temp*cos(omega*i):
  if (abs(temp) < .01) then     # How does one sum covariances out to infinity?
    break:                      #  I don't know so I cut off the loop once the covariance
  fi:                           #  gets below .01.
  if (YW = "" and i+xt[1] = time) then     # It also has to cut off if finding the covariance with Y[-1]
    break:
  fi:
od:
f := f/Pi:

RETURN(TSSubstitute(TS, f)):
end:


# Calls an assortment of procedures for exploratory time series analysis
ETSA := proc(TS :: list)
local xt, f, temp:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

print(`Time series model:`):
temp := TSSubstitute(TS, ARMAModel(xt[1], xt[2], 0)):
temp := Y[t] = temp:
print(temp):

print(`Autocorrelations:`):
print(TSPlot(TS, "YW")):

print(`Partial Autocorrelations:`):
print(TSPlot(TS, "partial", "YW")):

print(`Unit Roots Analysis:`):
UnitRoots(TS):

print(`Spectral Density Function:`):
f := SDF(TS, "YW"):
print(f):
print(plot(f,omega=0..Pi)):

print(`First Realization:`):
TSRealization(TS):
print(`Second Realization:`):
TSRealization(TS):

RETURN():
end:

# Generates and solves the Yule-Walker equations then finds the desired correlation or covariance
YuleWalker := proc(TS :: list, lag :: nonnegint)
local xt, phi, i, j, X, Y, temp, M, N, YW, variance, eq:
global rho:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (xt[4] = "time-dependent") then
  print(`ERROR(Yule-Walker): Error terms must be i.i.d.`):
  RETURN():
fi:

variance := Variance(xt[3]):
rho := array(0..xt[1]):
Y := {rho[0]}:
if (type(TS[1], list(realcons))) then
  M := convert([TS[1]], Matrix):
  phi := TS[1]:
else
  phi := array(1..xt[1]):
  for i from 1 to xt[1] do
    if (i = 1) then
      M := [phi[1]]:
      next:
    fi:
    M := [op(M), phi[i]]:
  od:
  M := convert(M, Matrix):
fi:
if (xt[2] = 0) then
  for i from 0 to xt[1] do
    for j from 1 to xt[1] do
      if (j = 1) then
	N := rho[abs(i-1)]:
	next:
      fi:
      N := N, rho[abs(i-j)]:
    od:
    N := <<N>>:
    if (member("var", [args]) and i = 0) then
      X := {rho[0] = (M.N)[1,1] + variance}:
      next:
    elif (i = 0) then
      X := {rho[0] = 1}:
      next:
    fi:
    X := {op(X), rho[i] = (M.N)[1,1]}:
    Y := {op(Y), rho[i]}:
  od:
else
  eq := YWEqConstruct(TS, 0, variance):
  X := {rho[0] = eq}:
  for i from 1 to xt[1] do
    eq := YWEqConstruct(TS, i, variance):
    X := {op(X), rho[i] = eq}:
    Y := {op(Y), rho[i]}:
  od:
fi:

temp := convert(solve(X, Y), 'string'):
temp := Split(temp[2..length(temp)-1], ","):
X := []:
for i from 1 to nops(temp) do
  X := [op(X), parse(temp[i])]:
od:
X := sort(X, (a, b) -> op(lhs(a)) < op(lhs(b))):
Y := []:
for i from 1 to nops(X) do
  Y := [op(Y), rhs(X[i])]:
od:

if (xt[2] > 0 and not member("var", [args])) then  # Transforms covariances to correlations in this special case
  for i from 2 to nops(Y) do
    Y[i] := Y[i] / Y[1]:
  od:
  Y[1] := 1:
fi:

if (lag <= xt[1]) then
  RETURN(Y[lag+1]):
else
  for i from xt[1]+1 to lag do
    temp := 0:
    for j from 1 to xt[1] do
      temp := temp + M[1,j]*Y[i-j+1]:
    od:
    Y := [op(Y), temp]:
  od:
  RETURN(Y[lag+1]):
fi:
end:

# To construct a system of equations for an ARMA(p, q), q>0 model, recursive substitution must be used
#  This is done here then called in the above procedure. See the section 2.4 of the associated paper for
#  more information. yecov = xi
YWEqConstruct := proc(TS :: list, lag :: nonnegint, variance)
local i, j, eq, phi, theta, yecov, xt:
global rho:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

rho := array(0..max(xt[1], xt[2])):
if (type(TS[1], list(realcons))) then
  phi := TS[1]:
else
  phi := array(1..xt[1]):
fi:
if (type(TS[2], list(realcons))) then
  theta := array(0..xt[2]):
  theta[0] := 1:
  for i from 1 to xt[2] do
    theta[i] := TS[2][i]:
  od:
else
  theta := array(0..xt[2]):
  theta[0] := 1:
fi:
eq := 0:
for i from 1 to xt[1] do
  eq := eq + phi[i]*rho[abs(i-lag)]:
od:

yecov := array(0..xt[2]):
yecov[0] := variance:
for i from 0 to xt[2] do
  if (i >= lag) then
    if (type(yecov[i-lag], indexed)) then
      yecov[i-lag] := theta[i-lag]*variance:
      for j from 1 to min(i-lag, xt[1]) do
	yecov[i-lag] := yecov[i-lag] + phi[j]*yecov[i-lag-j]:
      od:
    fi:
    eq := eq + theta[i]*yecov[i-lag]:
  fi:
od:

RETURN(eq):
end:

# Forecast an AR(p) model
TSForecast := proc(TS :: list, s :: posint, obs :: list(realcons))
local xt, epsilon, i, j, k, mean, var, temp, E, V, ystring, ylist, termlist, termmean, termvar, points, ub, lb, uconf, lconf:
global Y:

if (nops(TS) < 4) then
  xt := TSParse(TS):
else
  xt := TS:
fi:

if (xt[2] > 0) then
  print(`Forecasting is currently supported for AR(p) models only.`):
  RETURN():
fi:
if (xt[4] = "time-dependent") then
  print(`Forecasting is currently not supported for time-dependent models.`):
  RETURN():
fi:

Y := array(1..s):
epsilon := array(0..s):
for i from 0 to s-1 do
  Y[s-i] := epsilon[s-i]:
  for j from 1 to xt[1] do
    if (s-i-j < 1) then
      Y[s-i] := Y[s-i] + TS[1][j]*obs[abs(s-i-j-1)]:
    else
      Y[s-i] := Y[s-i] + TS[1][j]*Y[s-i-j]:
    fi:
  od:
od:

if (nargs > 3 and type(args[4], list)) then
  uconf := args[4][1]:
  lconf := args[4][2]:
else
  uconf := .95:
  lconf := .05:
  print(`Default confidence interval is 10%.`):
fi:
points := [seq([i, obs[nops(obs)-i+1]], i=1..nops(obs))]:

E := Mean(xt[3]):
V := Variance(xt[3]):
for k from 1 to s do
  mean := 0:
  var := 0:
  ystring := convert(Y[k], 'string'):
  ystring := SubstituteAll(ystring, "-", "+(-1)*"):
  if (ystring[1] = "+") then
    ystring := Insert(ystring, 0, "0"):
  fi:
  ylist := Split(ystring, "+"):
  for i from 1 to nops(ylist) do
    termmean := 1:
    termvar := 1:
    termlist := Split(ylist[i], "*"):
    for j from 1 to nops(termlist) do
      if (IsPrefix("epsilon", termlist[j])) then
	termmean := termmean * E:
	termvar := termvar * V:
      else
	termmean := termmean * parse(termlist[j]):
	termvar := termvar * parse(termlist[j])^2:
      fi:
    od:
    mean := mean + termmean:
    var := var + termvar:
  od:
  temp := NormalRV(mean, sqrt(var)):
  ub := IDF(temp, uconf):
  lb := IDF(temp, lconf):
  points := [op(points), [nops(obs)+k, mean], [nops(obs)+k, ub], [nops(obs)+k, lb]]:
od:

print(`Assuming error terms are normal for this forecast.`):
print(pointplot(points, axes=FRAME, labels=["t", "Y[t]"], xtickmarks=[seq(i, i=0..nops(obs)+s)])):
print(`Expected value and confidence interval are:`):
RETURN([mean, lb, ub]):
end: