# HW 3 Q 7 lag1-knn.03.R JLCY, 19 Dec 2013

############################################################ 

### This function implements K-NN lag-1 resampling It lags the input data
### vector by 1 into vectors Xt1- the previous time and Xt the current time
### neighbors are found on Xt1 and simulated from Xt

# length of the annual flow
N = length(annflow)

## Error: 找不到物件 'annflow'

N1 = N - 1

## Error: 找不到物件 'N'

Xt = annflow[2:N]

## Error: 找不到物件 'annflow'

Xt1 = annflow[1:N1]  # current state

## Error: 找不到物件 'annflow'

## number of simulations desired is nsim each simulation will be of length N

nsim = nsims  # you can change this

## Error: 找不到物件 'nsims'

simflow = matrix(0, ncol = nsim, nrow = N)

## Error: 找不到物件 'N'

simpdf = matrix(0, nrow = nsim, ncol = 100)

## Error: 找不到物件 'nsim'

## the weight metric to do the K-NN resampling
K = round(sqrt(N1))

## Error: 找不到物件 'N1'

W = 1:K

## Error: 找不到物件 'K'

W = 1/W

## Error: 找不到物件 'W'

W = W/sum(W)

## Error: 找不到物件 'W'

W = cumsum(W)

## Error: 找不到物件 'W'

for (i in 1:nsim) {

    # pick a flow at random to start

    if (Xt1[1] == 0) {
        if (Xt1[2] == 0) {
            simflow[1, i] = sample(pair00[, 1], 1)
        } else {
            simflow[1, i] = sample(pair01[, 1], 1)
        }
    } else {
        if (Xt1[2] == 0) {
            simflow[1, i] = sample(pair10[, 1], 1)
        } else {
            simflow[1, i] = sample(pair11[, 1], 1)
        }
    }

    xprev = simflow[1, i]  # xt1

    for (j in 1:N1) {

        # find the dist of xprev to all the N-1 years in the historical record order
        # the distances
        if (Xt1[j] == 0) {
            if (Xt[j] == 0) {
                Xt2 = pair00[, 2]
            } else {
                Xt2 = pair01[, 2]
            }
        } else {
            if (Xt[j] == 0) {
                Xt2 = pair10[, 2]
            } else {
                Xt2 = pair11[, 2]
            }

        }

        xdist = order(abs(xprev - Xt2))

        # select one of the nearest neighbor with the weight function above.

        xx = runif(1, 0, 1)
        xy = c(xx, W)
        xx = rank(xy)
        index = xdist[xx[1]]

        # the successor value of sampled index, which is in vector Xt
        simflow[j + 1, i] = Xt[index]
        xprev = Xt[index]
    }

    armean[i] = mean(simflow[, i])
    armax[i] = max(simflow[, i])
    armin[i] = min(simflow[, i])
    arvar[i] = var(simflow[, i])
    arskw[i] = skew(simflow[, i])

    # lag-1 correlation
    arcor[i] = cor(simflow[1:99, i], simflow[2:100, i])

    # annsim = apply(simflow[,i],1,sum)

    # simpdf[i,]=sm.density(simflow[,i],eval.points=100,display='none')$estimate
    simpdf[i, ] = sm.density(simflow[, i], eval.points = annxeval, display = "none")$estimate


}

## Error: 找不到物件 'nsim'