#install.packages("astsa")
require(astsa)
#help(package = astsa) 

eq = rowMeans(eqexp[,1:8]) # known earthquakes
ex = rowMeans(eqexp[,9:16]) # known explosions
event = eqexp[,17] # an event in Novaya Zenba (Russia); to be classified 

r = function(x) max(x) - min(x)  # The range of values in x
eq = eq/r(eq)
ex = ex/r(ex)
event = event/r(event)  # Scale to have range 1

par(mfrow=c(3,1))
ts.plot(eq, main = "Earthquake")
ts.plot(ex, main = "Explosion")
ts.plot(event, main = "Event")

# Each of these series is really two series - an early ("P") series of length 1024 and 
# a much later ("S) series of the same length.  The later signal travels through a different part of the earth's crust.
# The structure of the two types of signals tends to be different for earthquakes and explosions; they are compared 
# in order to differentiate between earthquakes and explosions.

# Split each of the series into its P and S components:

eqP = eq[1:1024]
eqS = eq[1025:2048]

exP = ex[1:1024]
exS = ex[1025:2048]

evP = event[1:1024]
evS = event[1025:2048]

win.graph()
par(mfrow=c(3,2))
ts.plot(eqP, main = "eqP")
ts.plot(eqS, main = "eqS")

ts.plot(exP, main = "exP")
ts.plot(exS, main = "exS")

ts.plot(evP, main = "evP")
ts.plot(evS, main = "evS")


# Can we use these data to decide whether the event was an earthquake or an expolosion?
# Look first at the spectra:

win.graph()
L = 21 # Chosen by trial and error
par(mfrow=c(3,2))
spec.eqP = spec.pgram(eqP,  kernel("daniell", (L-1)/2))
spec.eqS = spec.pgram(eqS,  kernel("daniell", (L-1)/2))
spec.exP = spec.pgram(exP,  kernel("daniell", (L-1)/2))
spec.exS = spec.pgram(exS,  kernel("daniell", (L-1)/2))
spec.evP = spec.pgram(evP,  kernel("daniell", (L-1)/2))
spec.evS = spec.pgram(evS,  kernel("daniell", (L-1)/2))

## Look at the coherences between "eq" and "ev", and between "ex" and "ev", for both P and S:

win.graph()
par(mfrow=c(2,2))
alpha = .05

spec.evP.eqP = spec.pgram(cbind(evP,eqP), kernel("daniell", (L-1)/2), plot.type = "coh", 
 main = "coherence; evP and eqP")
df = spec.evP.eqP$df
f = qf(1-alpha, 2, df-2) 
c = f/(.5*df-1+f)     
abline(h = c) 

spec.evS.eqS = spec.pgram(cbind(evS,eqS), kernel("daniell", (L-1)/2), plot.type = "coh", 
 main = "coherence; evS and eqS")
df = spec.evS.eqS$df
f = qf(1-alpha, 2, df-2) 
c = f/(.5*df-1+f)     
abline(h = c)

spec.evP.exP = spec.pgram(cbind(evP,exP), kernel("daniell", (L-1)/2), plot.type = "coh", 
 main = "coherence; evP and exP")
df = spec.evP.exP$df
f = qf(1-alpha, 2, df-2) 
c = f/(.5*df-1+f)     
abline(h = c)

spec.evS.exS = spec.pgram(cbind(evS,exS), kernel("daniell", (L-1)/2), plot.type = "coh", 
 main = "coherence; evS and exS")
df = spec.evS.exS$df
f = qf(1-alpha, 2, df-2) 
c = f/(.5*df-1+f)     
abline(h = c)


# The ranges of the S and P spectra, and in particular the ratios of these ranges, 
# are thought to be reasonable indicators

r.eq = c(r(spec.eqS$spec),r(spec.eqP$spec),r(spec.eqS$spec)/r(spec.eqP$spec))
r.ex = c(r(spec.exS$spec),r(spec.exP$spec),r(spec.exS$spec)/r(spec.exP$spec))
r.ev = c(r(spec.evS$spec),r(spec.evP$spec),r(spec.evS$spec)/r(spec.evP$spec))
rangemat = cbind(r.eq, r.ex, r.ev)
colnames(rangemat) = c("Eq", "Ex", "Ev")
rownames(rangemat) = c("S range", "P range", "ratio")
print(round(rangemat,4))

# Can we discriminate on the basis of impulse-response models predicting ev  
# from ex or eq, broken down by P and S?

M = 64

win.graph()
out.exP = LagReg(input = exP, output = evP, L, M, threshold=.05) # MSE =   0.01742299
win.graph()
out.exS = LagReg(input = exS, output = evS, L, M, threshold=.05) # MSE =   0.01985787


win.graph()
out.eqP = LagReg(input = eqP, output = evP, L, M, threshold=.05) # MSE =  0.01662466
win.graph()
out.eqS = LagReg(input = eqS, output = evS, L, M, threshold=.05) # MSE =  0.01936941


# Try filtering them first:
out.exP = SigExtract(exP, L, M, max.freq=.1)
out.eqP = SigExtract(eqP,  L, M, max.freq=.1)
out.evP = SigExtract(evP,  L, M, max.freq=.1)

out.exS = SigExtract(exS,  L, M, max.freq=.1)
out.eqS = SigExtract(eqS,  L, M, max.freq=.1)
out.evS = SigExtract(evS,  L, M, max.freq=.1)

# Now do LagReg again:

win.graph()
out2.exP = LagReg(input = na.omit(out.exP), output = na.omit(out.evP), L, M, threshold=.05) # MSE = 0.002264346
win.graph()
out2.exS = LagReg(input = na.omit(out.exS), output = na.omit(out.evS), L, M, threshold=.05) # MSE =   0.003665618

win.graph()
out2.eqP = LagReg(input = na.omit(out.eqP), output = na.omit(out.evP), L, M, threshold=.05) # MSE = 0.002199309
win.graph()
out2.eqS = LagReg(input = na.omit(out.eqS), output = na.omit(out.evS), L, M, threshold=.05) # MSE = 0.004097862


#graphics.off()